Skip to content

Commit bc8d640

Browse files
CopilotSierd
Copilot
and
Sierd
committed
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. Co-authored-by: Sierd <[email protected]>
1 parent 2031e2d commit bc8d640

File tree

2 files changed

+167
-2
lines changed

2 files changed

+167
-2
lines changed

‎aeolis/tests/test_wind.py‎

Lines changed: 165 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,165 @@
1+
"""
2+
This file is part of AeoLiS test suite. It uses the Pytest framework.
3+
4+
Tests for wind module functionality, specifically for wind shear velocity components.
5+
"""
6+
7+
import numpy as np
8+
import pytest
9+
10+
11+
class TestWindShearVelocityComponents:
12+
"""Test wind shear velocity component calculations"""
13+
14+
def test_ustars0_ustarn0_components_pure_s_direction(self):
15+
"""Test that ustars0 and ustarn0 correctly represent directional components.
16+
17+
For wind from 270 degrees (nautical, wind from west) with alfa=0:
18+
- Wind should be in pure s-direction (eastward)
19+
- ustars should be positive
20+
- ustarn should be 0
21+
- ustars0 should equal ustars (not the magnitude ustar)
22+
- ustarn0 should equal ustarn (which is 0, not the magnitude ustar)
23+
"""
24+
from aeolis.wind import interpolate, initialize
25+
from aeolis.model import ModelState
26+
27+
# Create minimal model state
28+
s = ModelState()
29+
30+
# Grid setup - simple 2x2 grid
31+
nx, ny = 3, 3
32+
s['x'] = np.zeros((ny, nx))
33+
s['y'] = np.zeros((ny, nx))
34+
s['zb'] = np.zeros((ny, nx))
35+
s['uw'] = np.zeros((ny, nx))
36+
s['udir'] = np.zeros((ny, nx))
37+
s['uws'] = np.zeros((ny, nx))
38+
s['uwn'] = np.zeros((ny, nx))
39+
s['ustars'] = np.zeros((ny, nx))
40+
s['ustarn'] = np.zeros((ny, nx))
41+
s['ustar'] = np.zeros((ny, nx))
42+
s['ustars0'] = np.zeros((ny, nx))
43+
s['ustarn0'] = np.zeros((ny, nx))
44+
s['ustar0'] = np.zeros((ny, nx))
45+
s['tau'] = np.zeros((ny, nx))
46+
s['taus'] = np.zeros((ny, nx))
47+
s['taun'] = np.zeros((ny, nx))
48+
s['tau0'] = np.zeros((ny, nx))
49+
s['taus0'] = np.zeros((ny, nx))
50+
s['taun0'] = np.zeros((ny, nx))
51+
s['hveg'] = np.zeros((ny, nx))
52+
53+
# Parameters
54+
p = {}
55+
p['alfa'] = 0.0 # Grid orientation (0 = aligned with north)
56+
p['kappa'] = 0.4 # von Karman constant
57+
p['z'] = 10.0 # Reference height (m)
58+
p['k'] = 0.001 # Roughness length (m)
59+
p['method_roughness'] = 'constant'
60+
p['rhoa'] = 1.25 # Air density (kg/m³)
61+
p['process_wind'] = False # We'll set wind manually
62+
p['wind_file'] = None
63+
64+
# Set constant wind from 270 degrees (from west, blowing east)
65+
# In nautical convention: 270° means wind FROM the west
66+
s['uw'][:, :] = 15.0 # 15 m/s
67+
s['udir'][:, :] = 270.0 # degrees nautical
68+
69+
# Run wind interpolation (which calculates ustars, ustarn, etc.)
70+
s = interpolate(s, p, t=0.0)
71+
72+
# For wind from 270° with alfa=0:
73+
# uws = -uw * sin((-alfa + udir) * pi/180) = -15 * sin(270°) = -15 * (-1) = 15
74+
# uwn = -uw * cos((-alfa + udir) * pi/180) = -15 * cos(270°) = -15 * 0 = 0
75+
76+
# Therefore:
77+
# ustars should be positive (proportional to uws)
78+
# ustarn should be 0 (proportional to uwn)
79+
80+
# Check that ustars is positive and ustarn is approximately 0
81+
assert np.all(s['ustars'] > 0), "ustars should be positive for wind from 270°"
82+
assert np.allclose(s['ustarn'], 0, atol=1e-10), "ustarn should be 0 for wind from 270°"
83+
84+
# The key test: ustars0 should equal ustars, not ustar
85+
# and ustarn0 should equal ustarn (0), not ustar
86+
assert np.allclose(s['ustars0'], s['ustars'], rtol=1e-10), \
87+
"ustars0 should equal ustars, not ustar magnitude"
88+
assert np.allclose(s['ustarn0'], s['ustarn'], rtol=1e-10), \
89+
"ustarn0 should equal ustarn (0), not ustar magnitude"
90+
91+
# Verify that ustar (magnitude) is not zero
92+
assert np.all(s['ustar'] > 0), "ustar magnitude should be positive"
93+
94+
# Verify that ustarn0 is NOT equal to ustar (the bug we're fixing)
95+
assert not np.allclose(s['ustarn0'], s['ustar'], rtol=1e-10), \
96+
"ustarn0 should NOT equal ustar magnitude (this was the bug)"
97+
98+
def test_ustars0_ustarn0_components_pure_n_direction(self):
99+
"""Test ustars0 and ustarn0 for wind purely in n-direction.
100+
101+
For wind from 0/360 degrees (nautical, wind from north) with alfa=0:
102+
- Wind should be in pure n-direction (southward)
103+
- ustars should be 0
104+
- ustarn should be positive
105+
"""
106+
from aeolis.wind import interpolate
107+
from aeolis.model import ModelState
108+
109+
# Create minimal model state
110+
s = ModelState()
111+
112+
# Grid setup
113+
nx, ny = 3, 3
114+
s['x'] = np.zeros((ny, nx))
115+
s['y'] = np.zeros((ny, nx))
116+
s['zb'] = np.zeros((ny, nx))
117+
s['uw'] = np.zeros((ny, nx))
118+
s['udir'] = np.zeros((ny, nx))
119+
s['uws'] = np.zeros((ny, nx))
120+
s['uwn'] = np.zeros((ny, nx))
121+
s['ustars'] = np.zeros((ny, nx))
122+
s['ustarn'] = np.zeros((ny, nx))
123+
s['ustar'] = np.zeros((ny, nx))
124+
s['ustars0'] = np.zeros((ny, nx))
125+
s['ustarn0'] = np.zeros((ny, nx))
126+
s['ustar0'] = np.zeros((ny, nx))
127+
s['tau'] = np.zeros((ny, nx))
128+
s['taus'] = np.zeros((ny, nx))
129+
s['taun'] = np.zeros((ny, nx))
130+
s['tau0'] = np.zeros((ny, nx))
131+
s['taus0'] = np.zeros((ny, nx))
132+
s['taun0'] = np.zeros((ny, nx))
133+
s['hveg'] = np.zeros((ny, nx))
134+
135+
# Parameters
136+
p = {}
137+
p['alfa'] = 0.0
138+
p['kappa'] = 0.4
139+
p['z'] = 10.0
140+
p['k'] = 0.001
141+
p['method_roughness'] = 'constant'
142+
p['rhoa'] = 1.25
143+
p['process_wind'] = False
144+
p['wind_file'] = None
145+
146+
# Set wind from 0/360 degrees (from north, blowing south)
147+
s['uw'][:, :] = 15.0
148+
s['udir'][:, :] = 0.0 # or 360.0
149+
150+
# Run wind interpolation
151+
s = interpolate(s, p, t=0.0)
152+
153+
# For wind from 0° with alfa=0:
154+
# uws = -uw * sin(0) = 0
155+
# uwn = -uw * cos(0) = -15
156+
157+
# Check components
158+
assert np.allclose(s['ustars'], 0, atol=1e-10), "ustars should be 0 for wind from 0°"
159+
assert np.all(s['ustarn'] < 0), "ustarn should be negative for wind from 0°"
160+
161+
# Check that the 0 components are preserved
162+
assert np.allclose(s['ustars0'], s['ustars'], rtol=1e-10), \
163+
"ustars0 should equal ustars (0)"
164+
assert np.allclose(s['ustarn0'], s['ustarn'], rtol=1e-10), \
165+
"ustarn0 should equal ustarn"

‎aeolis/wind.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -148,8 +148,8 @@ def interpolate(s, p, t):
148148
s = velocity_stress(s,p)
149149

150150
s['ustar0'] = s['ustar'].copy()
151-
s['ustars0'] = s['ustar'].copy()
152-
s['ustarn0'] = s['ustar'].copy()
151+
s['ustars0'] = s['ustars'].copy()
152+
s['ustarn0'] = s['ustarn'].copy()
153153

154154
s['tau0'] = s['tau'].copy()
155155
s['taus0'] = s['taus'].copy()

0 commit comments

Comments
 (0)