Allow user to specify universal Kolmogorov constant for TKE dissipation rate function (#406)

It was pointed out to me that a universal constant of 0.5 was hardcoded
into dolfyn when calculating dissipation rate using the Lumley and
Terray method. According to turbulence theory and historical
experimental validation, this is the value to use for the streamwise
velocity spectra. According to the same, the constant is different for
velocity spectra that are orthogonal to the streamwise direction -
usually taken to be somewhere between 0.65 and 0.69.

This PR allows the user to set this constant via function input.

Pope, Stephen B. "Turbulent flows." Measurement Science and Technology
12.11 (2001): 2020-2021. pp 231-232.

Author: jmcvey3

examples/data/dolfyn/test_data/Sig1000_tidal_bin.nc

@@ -835,19 +835,26 @@ def check_turbulence_cascade_slope(self, psd, freq_range=[0.2, 0.4]):
 
         return m, b
 
-    def dissipation_rate_LT83(self, psd, U_mag, freq_range=[0.2, 0.4], noise=None):
+    def dissipation_rate_LT83(
+        self, psd, U_mag, freq_range=[0.2, 0.4], k_constant=0.67, noise=None
+    ):
         """
         Calculate the TKE dissipation rate from the velocity spectra.
 
         Parameters
         ----------
         psd : xarray.DataArray (time,f)
-          The power spectral density from a single depth bin (range)
+          The power spectral density from the vertical beam and depth
+          bin (range)
         U_mag : xarray.DataArray (time)
-          The bin-averaged horizontal velocity (a.k.a. speed) from a single depth bin (range)
+          The bin-averaged horizontal velocity (a.k.a. speed) from a single
+          depth bin (range)
         f_range : iterable(2)
           The range over which to integrate/average the spectrum, in units
           of the psd frequency vector (Hz or rad/s)
+        k_constant : float or iterable(3)
+          Kolmogorov Constant (\\alpha in Notes section below) to use. Default
+          \\alpha is 0.67.
         noise : float or array-like
           Instrument noise level in same units as velocity. Typically
           found from `adp.turbulence.doppler_noise_level`.
@@ -864,10 +871,9 @@ def dissipation_rate_LT83(self, psd, U_mag, freq_range=[0.2, 0.4], noise=None):
 
         .. math:: S(k) = \\alpha \\epsilon^{2/3} k^{-5/3} + N
 
-        where :math:`\\alpha = 0.5` (1.5 for all three velocity
-        components), `k` is wavenumber, `S(k)` is the turbulent
-        kinetic energy spectrum, and `N' is the doppler noise level
-        associated with the TKE spectrum.
+        where :math:`\\alpha is the Kolmogorov constant (0.67 for vertical direction),
+        `k` is wavenumber, `S(k)` is the turbulent kinetic energy spectrum, and
+        `N' is the doppler noise level associated with the TKE spectrum.
 
         With :math:`k \\rightarrow \\omega / U`, then -- to preserve variance --
         :math:`S(k) = U S(\\omega)`, and so this becomes:
@@ -882,15 +888,19 @@ def dissipation_rate_LT83(self, psd, U_mag, freq_range=[0.2, 0.4], noise=None):
         by a random wave field". JPO, 1983, vol13, pp2000-2007.
         """
 
+        if not isinstance(psd, xr.DataArray):
+            raise TypeError("`psd` must be an instance of `xarray.DataArray`.")
         if len(psd.shape) != 2:
-            raise Exception("PSD should be 2-dimensional (time, frequency)")
+            raise Exception("`psd` should be 2-dimensional (time, frequency)")
         if len(U_mag.shape) != 1:
             raise Exception("U_mag should be 1-dimensional (time)")
         if not hasattr(freq_range, "__iter__") or len(freq_range) != 2:
             raise ValueError("`freq_range` must be an iterable of length 2.")
+        if np.size(k_constant) != 1:
+            raise ValueError("`k_constant` should be a single value.")
         if noise is not None:
             if np.shape(noise)[0] != np.shape(psd)[0]:
-                raise Exception("Noise should have same first dimension as PSD")
+                raise Exception("Noise should have same first dimension as `psd`")
         else:
             noise = np.array(0)
 
@@ -904,12 +914,15 @@ def dissipation_rate_LT83(self, psd, U_mag, freq_range=[0.2, 0.4], noise=None):
         idx = np.where((freq_range[0] < freq) & (freq < freq_range[1]))
         idx = idx[0]
 
+        # Set the correct magnitude whether the frequency is in Hz or rad/s
         if freq.units == "Hz":
             U = U_mag / (2 * np.pi)
         else:
             U = U_mag
 
-        a = 0.5
+        # Use the transverse value derived from the Kolmogorov constant
+        a = k_constant
+        # Calculate dissipation
         out = (psd[:, idx] * freq[idx] ** (5 / 3) / a).mean(axis=-1) ** (
             3 / 2
         ) / U.values

examples/data/dolfyn/test_data/vector_data01_bin.nc

@@ -339,9 +339,16 @@ def check_turbulence_cascade_slope(self, psd, freq_range=[6.28, 12.57]):
 
         return m, b
 
-    def dissipation_rate_LT83(self, psd, U_mag, freq_range=[6.28, 12.57], noise=None):
+    def dissipation_rate_LT83(
+        self,
+        psd,
+        U_mag,
+        freq_range=[6.28, 12.57],
+        k_constant=[0.5, 0.67, 0.67],
+        noise=None,
+    ):
         """
-        Calculate the dissipation rate from the PSD
+        Calculate the dissipation rate from the power spectral density of velocity.
 
         Parameters
         ----------
@@ -353,6 +360,12 @@ def dissipation_rate_LT83(self, psd, U_mag, freq_range=[6.28, 12.57], noise=None
           The range over which to integrate/average the spectrum, in units
           of the psd frequency vector (Hz or rad/s).
           Default = [6.28, 12.57] rad/s
+        k_constant : float or iterable(3)
+          Kolmogorov Constant (\\alpha in Notes section below) to use. If a
+          three dimensional PSD is provided, \\alpha defaults to [0.5, 0.67, 0.67];
+          i.e. 0.5 for the streamwise PSD and 0.67 for the transverse and vertical
+          PSDs. If the PSD is provided for a single velocity direction, \\alpha is
+          taken to be 0.5 unless otherwise specified.
         noise : float or array-like
           Instrument noise level in same units as velocity. Typically
           found from `adv.turbulence.calc_doppler_noise`.
@@ -369,10 +382,9 @@ def dissipation_rate_LT83(self, psd, U_mag, freq_range=[6.28, 12.57], noise=None
 
         .. math:: S(k) = \\alpha \\epsilon^{2/3} k^{-5/3} + N
 
-        where :math:`\\alpha = 0.5` (1.5 for all three velocity
-        components), `k` is wavenumber, `S(k)` is the turbulent
-        kinetic energy spectrum, and `N' is the doppler noise level
-        associated with the TKE spectrum.
+        where :math:`\\alpha is the Kolmogorov constant, `k` is wavenumber,
+        `S(k)` is the turbulent kinetic energy spectrum, and `N' is the
+        doppler noise level associated with the TKE spectrum.
 
         With :math:`k \\rightarrow \\omega / U`, then -- to preserve variance --
         :math:`S(k) = U S(\\omega)`, and so this becomes:
@@ -390,15 +402,19 @@ def dissipation_rate_LT83(self, psd, U_mag, freq_range=[6.28, 12.57], noise=None
         if not isinstance(psd, xr.DataArray):
             raise TypeError("`psd` must be an instance of `xarray.DataArray`.")
         if len(U_mag.shape) != 1:
-            raise Exception("U_mag should be 1-dimensional (time)")
+            raise Exception("U_mag should be 1-dimensional (time).")
         if len(psd["time"]) != len(U_mag["time"]):
-            raise Exception("`U_mag` should be from ensembled-averaged dataset")
+            raise Exception("`U_mag` should be from ensembled-averaged dataset.")
         if not hasattr(freq_range, "__iter__") or len(freq_range) != 2:
             raise ValueError("`freq_range` must be an iterable of length 2.")
-
+        # if the spectra are 1D, then the first dimension should be time (any length)
+        if (psd.shape[0] != 3) and (np.size(k_constant) != 1):
+            raise ValueError("`k_constant` should be a single value.")
+        elif (psd.shape[0] == 3) and (np.size(k_constant) != 3):
+            raise ValueError("`k_constant` should be an iterable of length 3.")
         if noise is not None:
-            if np.shape(noise)[0] != 3:
-                raise Exception("Noise should have same first dimension as velocity")
+            if np.shape(noise)[0] != np.shape(psd)[0]:
+                raise Exception("Noise should have same first dimension as `psd`.")
         else:
             noise = np.array([0, 0, 0])[:, None, None]
 
@@ -412,12 +428,20 @@ def dissipation_rate_LT83(self, psd, U_mag, freq_range=[6.28, 12.57], noise=None
         idx = np.where((freq_range[0] < freq) & (freq < freq_range[1]))
         idx = idx[0]
 
+        # Set the correct magnitude whether the frequency is in Hz or rad/s
         if freq.units == "Hz":
             U = U_mag / (2 * np.pi)
         else:
             U = U_mag
 
-        a = 0.5
+        # Set Kolmogorov constant
+        a = np.array(k_constant)
+        if psd.shape[0] == 3:
+            a = a[:, None, None]  # stack properly
+        else:
+            a = np.squeeze(k_constant)
+
+        # Calculate dissipation
         out = (psd.isel(freq=idx) * freq.isel(freq=idx) ** (5 / 3) / a).mean(
             axis=-1
         ) ** (3 / 2) / U

mhkit/dolfyn/adp/turbulence.py

@@ -143,7 +143,7 @@ def test_adcp_turbulence(self):
         dat.velds.rotate2("beam")
 
         tdat["psd"] = bnr.power_spectral_density(
-            dat["vel"].isel(dir=2, range=len(dat["range"]) // 2), freq_units="Hz"
+            dat["vel_b5"].isel(range_b5=len(dat["range_b5"]) // 2), freq_units="Hz"
         )
         tdat["noise"] = bnr.doppler_noise_level(tdat["psd"], pct_fN=0.8)
         tdat["stress_vec4"] = bnr.reynolds_stress_4beam(
@@ -179,11 +179,11 @@ def test_adcp_turbulence(self):
         ) = bnr.dissipation_rate_SF(dat["vel"].isel(dir=2), r_range=[1, 5])
 
         slope_check = bnr.check_turbulence_cascade_slope(
-            tdat["psd"].mean("time"), freq_range=[0.4, 4]
+            tdat["psd"].mean("time_b5"), freq_range=[0.4, 4]
         )
         # Check noise subtraction in psd function
         tdat["psd_noise"] = bnr.power_spectral_density(
-            dat["vel"].isel(dir=2, range=len(dat["range"]) // 2),
+            dat["vel_b5"].isel(range_b5=len(dat["range_b5"]) // 2),
             freq_units="Hz",
             noise=0.01,
         )