"""Physical turbulence, drawing ``sigma_u`` for each sampled case.
A turbulence model is a link of the site-condition chain. It draws ``sigma_u``
in ``m/s`` from a quantile given the case's wind speed, surface and stability,
and exposes its conditional distribution of turbulence intensity so wind
resources can build TI roses from the same model they sample. The coverage
band of Dimitrov et al. is part of :class:`oowl.resources.Dimitrov`.
References
----------
de Maré, M. and Mann, J. (2016), On the space-time structure of sheared
turbulence, Boundary-Layer Meteorol. 160, 453-474.
https://doi.org/10.1007/s10546-016-0143-z
Kelly, M. (2018), From standard wind measurements to spectral characterization,
Wind Energ. Sci. 3, 533-543. https://doi.org/10.5194/wes-3-533-2018
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import TYPE_CHECKING, Any
import numpy as np
from .._repr import _FieldRepr
from .._validation import _validate_real, _validate_real_array
from .stability import (
_VON_KARMAN,
_checked_positive_speed,
_checked_quantiles,
_profile,
)
if TYPE_CHECKING:
from collections.abc import Sequence
from numpy.typing import ArrayLike
from .stability import KellyGryning
# Midpoint nodes over the stability quantile, so averages over them average
# over the 1/L the climate draws.
_STABILITY_NODES = (np.arange(256) + 0.5) / 256
# Step of the log-scale grid of a distributed coefficient, as a fraction of its
# coefficient of variation.
_GRID_STEP = 0.01
# Entries of one block of the grid tables, about 16 MB.
_BLOCK_SIZE = 2**21
[docs]
@dataclass(frozen=True, repr=False)
class SurfaceLayer(_FieldRepr):
"""Surface-layer turbulence, ``sigma_u = C u_star`` with the diabatic ``u_star``.
The friction velocity comes from the Monin-Obukhov profile at the case's
own stability, ``u_star = kappa U / [ln(z / z0) - Psi(z / L) + Psi(z0 / L)]``,
so unstable air is more turbulent and stable air less. Beyond
:data:`~oowl.resources.profile.STABLE_LIMIT` the stable gradient is held at the
critical Richardson number, so ``Psi`` grows logarithmically and
turbulence keeps falling slowly into very stable air.
``coefficient`` is ``C``, one number or a frozen ``scipy.stats``
distribution with bounded non-negative support, drawn independently for
every case. Kelly (2018) found ``2.3`` on average over flat land at
Høvsøre from 4 to 25 m/s, the default. ``C = 2.5`` gives the neutral rule
of thumb ``TI = 1 / ln(z / z0)``.
``intensity`` optionally sets the level, one value or one per sector. A
resource then scales ``sigma_u`` in each sector so its mean turbulence
intensity over 4 to 25 m/s is that value, keeping the spread over speed and
stability. The level depends on the site's winds, so the resource applies
it and the methods here give the unscaled model.
"""
#: ``sigma_u / u_star``, a number or a bounded ``scipy.stats`` distribution.
coefficient: Any = 2.3
#: Mean turbulence intensity over 4 to 25 m/s, one or one per sector.
intensity: float | Sequence[float] | np.ndarray | None = None
def __post_init__(self):
if self.intensity is not None:
values = _validate_real_array("intensity", self.intensity)
if values.ndim > 1 or values.size == 0:
raise ValueError("intensity must be one value or one per sector")
if not np.all(np.isfinite(values)) or np.any(values <= 0.0):
raise ValueError("intensity must contain positive finite values")
object.__setattr__(
self,
"intensity",
float(values) if values.ndim == 0 else tuple(values.tolist()),
)
coefficient = self.coefficient
if hasattr(coefficient, "ppf"):
low, high = coefficient.support()
if not np.isfinite(low) or not np.isfinite(high) or low < 0.0:
raise ValueError(
"coefficient distribution must have bounded non-negative support"
)
return
if isinstance(coefficient, (str, bytes, bool, np.bool_)):
value = None
else:
try:
value = float(coefficient)
except (TypeError, ValueError):
value = None
if value is None:
raise TypeError(
"coefficient must be a positive number or a scipy.stats "
f"distribution, not {type(coefficient).__name__}"
)
if not np.isfinite(value) or value <= 0.0:
raise ValueError("coefficient must be positive and finite")
object.__setattr__(self, "coefficient", value)
def _coefficient_ppf(self, quantiles: np.ndarray) -> np.ndarray:
if hasattr(self.coefficient, "ppf"):
return np.asarray(self.coefficient.ppf(quantiles), dtype=float)
return np.full(quantiles.shape, self.coefficient)
def _coefficient_support(self) -> tuple[float, float]:
if hasattr(self.coefficient, "support"):
low, high = self.coefficient.support()
return float(low), float(high)
return self.coefficient, self.coefficient
def _coefficient_mean(self) -> float:
if hasattr(self.coefficient, "mean"):
return float(self.coefficient.mean())
return self.coefficient
[docs]
def friction_velocity(
self,
wind_speed: ArrayLike,
*,
height: float,
roughness_length: ArrayLike,
inverse_obukhov_length: ArrayLike,
) -> np.ndarray:
"""Return the diabatic ``u_star`` in ``m/s`` at each case."""
u = _checked_positive_speed(wind_speed)
denominator = _profile(height, roughness_length, inverse_obukhov_length)
return _VON_KARMAN * u / denominator
[docs]
def sample(
self,
quantiles: ArrayLike,
*,
wind_speed: ArrayLike,
height: float,
roughness_length: ArrayLike,
inverse_obukhov_length: ArrayLike,
) -> np.ndarray:
"""Return ``sigma_u`` in ``m/s``, ``C`` at each quantile times ``u_star``."""
u_star = self.friction_velocity(
wind_speed,
height=height,
roughness_length=roughness_length,
inverse_obukhov_length=inverse_obukhov_length,
)
q = _checked_quantiles(quantiles, u_star.shape)
return self._coefficient_ppf(q) * u_star
def _intensity_scale(
self, wind_speed, height, roughness_length, stability
) -> np.ndarray:
"""Return ``kappa / [ln(z/z0) - Psi]`` over the stability quantile nodes.
The result has one row per wind speed and one column per node, so
``TI = C * scale`` integrates over the stability climate exactly drawn.
"""
u = _checked_positive_speed(wind_speed)
roughness = np.broadcast_to(np.asarray(roughness_length, dtype=float), u.shape)
inverse_length = stability._quantile_table(
u, _STABILITY_NODES, height=height, roughness_length=roughness
)
return _VON_KARMAN / _profile(height, roughness[:, None], inverse_length)
[docs]
def intensity_range(
self,
wind_speed: ArrayLike,
*,
height: float,
roughness_length: ArrayLike,
stability: KellyGryning,
) -> tuple[np.ndarray, np.ndarray]:
"""Return the lowest and highest turbulence intensity at each speed.
The extremes pair the coefficient support with the most stable and the
most unstable ``1/L`` the climate can draw, its quantiles 0 and 1.
"""
u = _checked_positive_speed(wind_speed)
low, high = self._coefficient_support()
scales = [
_VON_KARMAN
/ _profile(
height,
roughness_length,
stability.sample(
u,
np.full(u.shape, quantile),
height=height,
roughness_length=roughness_length,
),
)
for quantile in (0.0, 1.0)
]
return low * np.minimum(*scales), high * np.maximum(*scales)
[docs]
def intensity_cdf(
self,
wind_speed: ArrayLike,
values: ArrayLike,
*,
height: float,
roughness_length: ArrayLike,
stability: KellyGryning,
) -> np.ndarray:
"""Return ``P(TI < value | U)``, one row per speed and column per value.
Given the stability, ``TI`` lies below a value when ``C`` lies below it
divided by the scale, so the coefficient CDF is averaged over the
climate, drawn once for every value. A fixed ``C`` counts exactly. A
distributed ``C`` interpolates its CDF over a grid in the log of the
scale, so it is evaluated a few hundred times per value rather than
at every speed and stability. The kinks of the CDF at the ends of the
support are summed exactly, so a piecewise-linear CDF such as the
uniform is exact to rounding. A truncated normal is within about
``2e-9`` and a density with a kink inside its support within about
``1e-7``, while a density unbounded at an end, such as a beta with a
shape below one, errs by up to about ``1e-4`` near it. The result is
clipped to ``[0, 1]`` and made non-decreasing in the value.
"""
scale = self._intensity_scale(wind_speed, height, roughness_length, stability)
values = np.asarray(values, dtype=float)
if not hasattr(self.coefficient, "cdf"):
# P(C < x) of a fixed C, so bins closed below count it.
return np.stack(
[np.mean(value / scale > self.coefficient, axis=1) for value in values],
axis=1,
)
return _distributed_intensity_cdf(self.coefficient, scale, values)
[docs]
def mean_intensity(
self,
wind_speed: ArrayLike,
*,
height: float,
roughness_length: ArrayLike,
stability: KellyGryning,
) -> np.ndarray:
"""Return the mean turbulence intensity at each wind speed."""
scale = self._intensity_scale(wind_speed, height, roughness_length, stability)
return self._coefficient_mean() * np.mean(scale, axis=1)
[docs]
@dataclass(frozen=True, repr=False)
class Mann(_FieldRepr):
"""Mann turbulence parameters from each case's turbulence and shear.
The length scale follows Kelly (2018), ``l_M = sigma_u / (dU/dz)``, with
``dU/dz = U alpha / z`` from the case's own profile at the reference
height, so it grows in unstable air and shrinks in stable air. With
:class:`SurfaceLayer` it is ``C kappa z / phi_m(z/L)``. The anisotropy is
``anisotropy``, about 3 in measurements (de Maré and Mann, 2016, and Kelly,
2018) against the 3.9 of IEC 61400-1. Pass it as a resource's ``mann``.
"""
#: Eddy-lifetime parameter ``Gamma``, from 0 isotropic to about 5.
anisotropy: float = 3.0
#: Sobol dimensions the link draws.
dimensions = 0
def __post_init__(self):
value = _validate_real("anisotropy", self.anisotropy)
if not np.isfinite(value) or not 0.0 <= value <= 5.0:
raise ValueError("anisotropy must be finite and in [0, 5]")
object.__setattr__(self, "anisotropy", value)
[docs]
def sample(
self,
quantiles: ArrayLike,
*,
wind_speed: ArrayLike,
sigma_u: ArrayLike,
shear_exponent: ArrayLike,
height: float,
) -> dict[str, np.ndarray]:
"""Return ``mann_length_scale`` in meters and ``mann_anisotropy``.
``quantiles`` holds no rows, since the link is deterministic.
"""
u = _checked_positive_speed(wind_speed)
gradient = u * np.asarray(shear_exponent, dtype=float) / float(height)
return {
"mann_length_scale": np.asarray(sigma_u, dtype=float) / gradient,
"mann_anisotropy": np.full(u.shape, self.anisotropy),
}
def _distributed_intensity_cdf(
coefficient, scale: np.ndarray, values: np.ndarray
) -> np.ndarray:
"""Return ``P(C < value / scale)`` averaged over each row of ``scale``.
The CDF ``F`` of a continuous ``C`` with support ``[lo, hi]`` splits into
``alpha max(c - lo, 0) - beta max(c - hi, 0)``, whose row means are sums
over the scales, and a remainder with a continuous slope. Each log scale
takes the remainder by cubic interpolation from a uniform grid, whose step
is :data:`_GRID_STEP` times the coefficient of variation of ``C``, so the
CDF is evaluated only where the grid meets the support. ``alpha`` and
``beta`` are the mean slopes of ``F`` over a short step inside each end,
which keeps them finite where the end density is not.
"""
n_rows = scale.shape[0]
result = np.full((n_rows, values.size), np.nan)
known = ~np.isnan(values)
low, high = (float(end) for end in coefficient.support())
unique, inverse = np.unique(values[known], return_inverse=True)
# Below every scale's low end the CDF is zero, above every high end one.
smallest = scale.min(axis=1, keepdims=True)
largest = scale.max(axis=1, keepdims=True)
table = (unique >= high * largest).astype(float)
partial = (unique > low * smallest) & ~(unique >= high * largest)
columns = np.flatnonzero(partial.any(axis=0))
if columns.size:
within = unique[columns]
estimate = _interpolated_cdf(coefficient, low, high, scale, within)
table[:, columns] = np.where(partial[:, columns], estimate, table[:, columns])
# The interpolation can stray by rounding past zero, one and the order of
# the values, which the true CDF keeps.
np.minimum(np.maximum.accumulate(np.maximum(table, 0.0), axis=1), 1.0, out=table)
result[:, known] = table[:, inverse]
return result
def _interpolated_cdf(coefficient, low, high, scale, values) -> np.ndarray:
"""Return the row means of ``F(value / scale)`` for sorted ``values``."""
n_rows, n_nodes = scale.shape
step = _GRID_STEP * float(coefficient.std()) / float(coefficient.mean())
width = (high - low) * step
lower_slope = float(coefficient.cdf(low + width)) / width
upper_slope = float(coefficient.sf(high - width)) / width
# Four-node cubic Lagrange weights of each log scale on the grid.
log_scale = np.log(scale)
origin = log_scale.min() - 1.5 * step
position = (log_scale - origin) / step
node = np.floor(position).astype(np.intp)
offset = position - node
below, across = offset * (offset - 1.0), (offset + 1.0) * (offset - 2.0)
weights = np.stack(
[
below * (2.0 - offset) / 6.0,
across * (offset - 1.0) / 2.0,
across * offset / -2.0,
below * (offset + 1.0) / 6.0,
],
axis=1,
)
n_grid = int(node.max()) + 3
grid = np.exp(-(origin + step * np.arange(n_grid)))
shifts = np.arange(-1, 3)[:, None]
estimate = np.empty((n_rows, values.size))
block = max(1, _BLOCK_SIZE // n_grid)
for first in range(0, values.size, block):
# The remainder at every grid node, the CDF taken in the support.
part = slice(first, first + block)
ratio = values[None, part] * grid[:, None]
remainder = (ratio >= high).astype(float)
inside = (ratio > low) & (ratio < high)
remainder[inside] = coefficient.cdf(ratio[inside])
remainder -= lower_slope * np.maximum(ratio - low, 0.0)
remainder += upper_slope * np.maximum(ratio - high, 0.0)
for start in range(0, n_rows, block):
stop = min(start + block, n_rows)
local = np.arange(stop - start)[:, None] * n_grid + node[start:stop]
interpolation = np.bincount(
(local[:, None, :] + shifts).ravel(),
weights=weights[start:stop].ravel(),
minlength=(stop - start) * n_grid,
).reshape(-1, n_grid)
estimate[start:stop, part] = interpolation @ remainder
# Sums of max(value / scale - end, 0) over each row, from the count and the
# sum of the inverse scales above every threshold end / value.
thresholds = np.concatenate([low / values, high / values])
order = np.argsort(thresholds)
inverse_scale = 1.0 / scale
rank = np.searchsorted(thresholds[order], inverse_scale)
flat = (np.arange(n_rows)[:, None] * (thresholds.size + 1) + rank).ravel()
above = []
for weight in (None, inverse_scale.ravel()):
histogram = np.bincount(
flat, weights=weight, minlength=n_rows * (thresholds.size + 1)
).reshape(n_rows, -1)
tail = np.cumsum(histogram[:, ::-1], axis=1)[:, ::-1][:, 1:]
above.append(np.empty_like(tail))
above[-1][:, order] = tail
count, total = above
for slope, end, part in (
(lower_slope, low, slice(0, values.size)),
(-upper_slope, high, slice(values.size, None)),
):
estimate += slope * (values * total[:, part] - end * count[:, part])
return estimate / n_nodes
def _resolve_turbulence(turbulence):
"""Return the turbulence model accepted by a wind resource."""
if turbulence is None:
return SurfaceLayer()
if isinstance(turbulence, SurfaceLayer):
return turbulence
raise TypeError("turbulence must be a SurfaceLayer model")
