Fix ustarn calculation: initialization and FFT shear formula bugs (#265)

[Sierd]

• Initial plan

• Fix ustars0 and ustarn0 initialization bug in wind.py

Fixed bug where ustars0 and ustarn0 were incorrectly set to ustar magnitude instead of their respective directional components ustars and ustarn.

• input files for debugging

• Fix missing division in dtauy FFT shear calculation

The dtauy_t formula in the FFT shear method was missing the division by sc_kv(0., 2.*sqrt(2.)*sigma), causing incorrect y-direction shear stress perturbations. This resulted in non-zero ustarn values even when the bed had no y-direction variability.

The formula now matches the structure of dtaux_t calculation.

• Fix frequency array alignment in FFT shear calculation

The kx and ky frequency arrays were misaligned with the FFT output. The code was creating frequency arrays with fftfreq(n+1)[1:] which skips the DC component, but the FFT includes all frequencies. This caused incorrect mapping between frequencies and FFT coefficients, leading to wrong shear perturbations.

Fixed by using fftfreq(n) which properly matches the FFT output dimensions and includes the DC component at the correct index.

• Handle division by zero at DC component in FFT shear

The previous fix included the DC component (kx=0, ky=0) in the frequency arrays, which caused division by zero errors in the shear calculations.

Fixed by:

• Using safe division with np.where to replace zeros with 1.0 temporarily
• Explicitly setting DC component of perturbations to 0 after calculation
• Applying same fix to filter_highfrequencies function

The DC component represents the mean value and doesn't contribute to perturbations, so setting it to zero is physically correct.

• Use masked computation to avoid invalid values in FFT shear

The previous fix still allowed invalid values when kx=0 (but ky≠0) because sigma depends on kx. When kx=0, sigma=0 which causes issues in Bessel function ratios.

New approach:

• Create a mask for valid computations: (k > 0) & (abs(kx) > 0)
• Initialize perturbation arrays with zeros
• Only compute perturbations for valid frequency combinations
• All invalid combinations (kx=0 or k=0) automatically get zero perturbation

This is physically correct: perturbations at kx=0 (pure y-direction frequencies) should be zero for the x-direction flow formulation.

• Simplify FFT shear edge case handling

Simplified the masked computation approach to a cleaner implementation:

• Use np.where for safe division (replace zeros with 1.0 temporarily)
• Compute formulas normally with safe arrays
• Apply invalid_mask at the end to zero out problematic regions

This achieves the same result with much simpler, more readable code.

• Update shear.py

reverted the suggested division in dtauy_t calculation

• deleted wind tests

• shear patch

---

[Sierd]

@bartvanwesten can you pull the latest main into this...