Draft of lowfreq model in package now

.github/workflows/lint.yml

@@ -7,8 +7,13 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v4
-      - uses: psf/black@stable
+
+      - uses: chartboost/ruff-action@v1
+        with:
+          src: "./drdmannturb"
+          args: "--fix"
+
+      - uses: jpetrucciani/mypy-check@master
         with:
-          # options: "--check --verbose"
           src: "./drdmannturb"
-          version: "~= 22.3.0"
+          args: "--ignore-missing-imports"

.gitignore

@@ -1,20 +1,20 @@
-# csv data 
+# csv data
 *.csv
 
-# model pickles 
-*.pkl 
+# model pickles
+*.pkl
 
 # overrides for data used in examples
-# necessary for interpolation data example 
-!/docs/source/data/*.csv 
-# necessary for model saving, loading, and viz examples 
+# necessary for interpolation data example
+!/docs/source/data/*.csv
+# necessary for model saving, loading, and viz examples
 !/docs/source/results/*.pkl
 
 docs/source/auto_examples
 
 examples/runs/*
 
-# runtimes for sphinx-gallery 
+# runtimes for sphinx-gallery
 docs/source/sg_execution_times.rst
 
 test/io/runs
@@ -71,6 +71,7 @@ coverage.xml
 *.py,cover
 .hypothesis/
 .pytest_cache/
+.ruff_cache/
 
 # Translations
 *.mo
@@ -152,6 +153,7 @@ dmypy.json
 # Others
 **.vscode
 **.DS_Store
+examples-unrendered/2d_plus_3d.ipynb
 
 ## Data/ Paraview Files
 
@@ -162,4 +164,4 @@ runs/
 
 ## Open Journals compilation intermediates and results
 paper/jats
-paper/paper.pdf
+paper/paper.pdf

drdmannturb/fluctuation_generation/lowfreq_syed_mann.py

@@ -0,0 +1,193 @@
+"""
+This module implements the Syed-Mann (2024) low-frequency wind fluctuation model.
+"""
+
+from typing import Optional
+
+import numpy as np
+from scipy import integrate
+
+
+def _compute_kappa(k1: float, k2: float, psi: float) -> float:
+    """
+    Subroutine to compute the horizontal wavevector :math:`\kappa`, defined by
+
+    .. math::
+        \kappa = \sqrt{2(k_1^2 \cos^2(\psi) + k_2^2 \sin^2(\psi))}
+
+    Parameters
+    ----------
+    k1 : float
+        Wavenumber k1
+
+    k2 : float
+        Wavenumber k2
+
+    psi : float
+        "Anisotropy parameter" angle :math:`\psi`, in radians
+
+    Returns
+    -------
+    float
+        Computed kappa value
+    """
+
+    return np.sqrt(2.0 * ((k1**2) * np.cos(psi) ** 2 + (k2**2) * np.sin(psi) ** 2))
+
+
+def _compute_E(kappa: float, c: float, L2D: float, z_i: float) -> float:
+    """
+    Subroutine to compute the energy spectrum :math:`E(\kappa)` with the attenuation factor,
+    defined by
+
+    .. math::
+        E(\kappa) = \frac{c \kappa^3}{(L_{2\textrm{D}}^{-2} + \kappa^2)^{7/3}} \cdot
+        \frac{1}{1 + \kappa^2 z_i^2}
+
+    Parameters
+    ----------
+    kappa : float
+        Replacement "wavenumber" :math:`\kappa`
+
+    c : float
+        Scaling factor :math:`c` used to correct the variance
+
+    L2D : float
+        Length scale :math:`L_{2\textrm{D}}`
+    """
+    if np.isclose(kappa, 0.0):
+        return 0.0
+
+    denom = (1.0 / (L2D**2) + kappa**2) ** (7.0 / 3.0)
+    atten = 1.0 / (1.0 + (kappa * z_i) ** 2)
+    return c * (kappa**3) / denom * atten
+
+
+def _estimate_c(sigma2: float, L2D: float, z_i: float) -> float:
+    """
+    Subroutine to estimate the scaling factor :math:`c` from the target variance :math:`\sigma^2`.
+
+    This is achieved by approximating the integral of :math:`E(\kappa)` from :math:`\kappa=0` to
+    :math:`\infty` by quadrature, since
+    .. math::
+        \int_0^\infty E(\kappa)
+        = c \int_0^\infty \frac{\kappa^3}{(L_{2\textrm{D}}^{-2} + \kappa^2)^{7/3}} \cdot
+            \frac{1}{1 + \kappa^2 z_i^2}
+        = \sigma^2
+
+    Parameters
+    ----------
+    sigma2 : float
+        Target variance :math:`\sigma^2`
+
+    L2D : float
+        Length scale :math:`L_{2\textrm{D}}`
+
+    z_i : float
+        Height :math:`z_i`
+    """
+
+    def integrand(kappa: float) -> float:
+        return kappa**3 / ((1.0 / (L2D**2) + kappa**2) ** (7.0 / 3.0)) * (1.0 / (1.0 + (kappa * z_i) ** 2))
+
+    val, err = integrate.quad(integrand, 0, np.inf)
+
+    return sigma2 / val
+
+
+def generate_2D_lowfreq(
+    Nx: int,
+    Ny: int,
+    L1: float,
+    L2: float,
+    psi_degs: float,
+    sigma2: float,
+    L2D: float,
+    z_i: float,
+    c: Optional[float] = None,
+) -> np.ndarray:
+    """
+    Generates the 2D low-frequency wind fluctuation component of the Syed-Mann (2024) 2D+3D model.
+
+    Parameters
+    ----------
+    Nx : int
+        Number of grid points in the x-direction
+    Ny : int
+        Number of grid points in the y-direction
+    L1 : float
+        Length of the domain in the x-direction
+    L2 : float
+        Length of the domain in the y-direction
+    psi_degs : float
+        "Anisotropy parameter" angle :math:`\psi`, in degrees
+    sigma2 : float
+        Target variance :math:`\sigma^2`
+    L2D : float
+        Length scale :math:`L_{2\textrm{D}}`
+    z_i : float
+        Height :math:`z_i`
+    c : float
+        Scaling factor :math:`c` to use for the energy spectrum. If not provided, it is
+        estimated by quadrature from the provided target variance :math:`\sigma^2`.
+
+    Returns
+    -------
+    np.ndarray
+        Generated 2D low-frequency wind fluctuation component. This is `Nx` by `Ny` by 2,
+        where the third dimension is the u- (longitudinal) and v-components (transverse).
+
+        TODO ^^
+    """
+
+    assert 0 < psi_degs and psi_degs < 90, "Anisotropy parameter psi_degs must be between 0 and 90 degrees"
+
+    psi = np.deg2rad(psi_degs)
+
+    if c is None:
+        c = _estimate_c(sigma2, L2D, z_i)
+
+    dx = L1 / Nx
+    dy = L2 / Ny
+
+    kx_arr = 2.0 * np.pi * np.fft.fftfreq(Nx, d=dx)
+    ky_arr = 2.0 * np.pi * np.fft.fftfreq(Ny, d=dy)
+    kx_arr = np.fft.fftshift(kx_arr)
+    ky_arr = np.fft.fftshift(ky_arr)
+
+    amp2 = np.zeros((Nx, Ny), dtype=np.float64)
+
+    factor_16 = (2.0 * np.pi**2) / L1
+
+    for ix in range(Nx):
+        for iy in range(Ny):
+            kx = kx_arr[ix]
+            ky = ky_arr[iy]
+
+            kappa = _compute_kappa(kx, ky, psi)
+            E_val = _compute_E(kappa, c, L2D, z_i)
+
+            if kappa < 1e-12:
+                phi_11 = 0.0
+            else:
+                phi_11 = E_val / (np.pi * kappa)
+
+            amp2[ix, iy] = factor_16 * phi_11
+
+    Uhat = np.zeros((Nx, Ny), dtype=np.complex128)
+
+    for ix in range(Nx):
+        for iy in range(Ny):
+            amp = np.sqrt(amp2[ix, iy])
+            phase = (np.random.normal() + 1j * np.random.normal()) / np.sqrt(2.0)
+            Uhat[ix, iy] = amp * phase
+
+    Uhat_unshift = np.fft.ifftshift(Uhat, axes=(0, 1))
+    u_field_complex = np.fft.ifft2(Uhat_unshift, s=(Nx, Ny))
+    u_field = np.real(u_field_complex)
+
+    var_now = np.var(u_field)
+    if var_now > 1e-12:
+        u_field *= np.sqrt(sigma2 / var_now)
+
+    return u_field

examples-unrendered/2d_plus_3d_prototype.py

@@ -0,0 +1,59 @@
+import numpy as np
+from low_freq_prototype import generate_2D_lowfreq_approx
+
+from drdmannturb.fluctuation_generation import GenerateFluctuationField
+
+# DRDMT 3D turbulence parameters.
+Type_Model = "Mann"
+
+z0 = 0.02
+zref = 90
+uref = 11.4
+ustar = uref * 0.41 / np.log(zref / z0)
+plexp = 0.2
+windprofiletype = "LOG"
+
+L_3D = 50
+Gamma = 2.5
+sigma = 0.01
+
+Lx = 60_000  # [m] = 60 km
+Ly = 15_000  # [m] = 15 km
+Lz = 5_000  # [m]
+
+grid_dimensions = np.array([Lx, Ly, Lz])
+grid_levels = np.array([6, 4, 4])
+
+nBlocks = 1
+
+seed = None
+
+# Generate 3d Mann box
+gen_mann = GenerateFluctuationField(
+    ustar,
+    zref,
+    grid_dimensions,
+    grid_levels,
+    length_scale=L_3D,
+    time_scale=Gamma,
+    energy_spectrum_scale=sigma,
+    model=Type_Model,
+    seed=seed,
+)
+
+fluctuation_field = gen_mann.generate_fluctuation_field()
+
+spacing = tuple(grid_dimensions / (2.0**grid_levels + 1))
+
+
+# 2D field parameters
+L_2D = 15_000.0
+sigma2 = 0.6
+z_i = 500.0
+psi_degs = 43.0
+
+L1, L2 = grid_dimensions[:2]
+Nx, Ny = 2 ** grid_levels[:2]
+
+_, _, u_field = generate_2D_lowfreq_approx(Nx, Ny, L1, L2, psi_degs, sigma2, L_2D, z_i)
+_, _, v_field = generate_2D_lowfreq_approx(Nx, Ny, L1, L2, psi_degs, sigma2, L_2D, z_i)

examples-unrendered/low_freq_prototype.py

@@ -7,7 +7,7 @@ def compute_kappa(kx, ky, psi):
     cos2 = np.cos(psi) ** 2
     sin2 = np.sin(psi) ** 2
 
-    return np.sqrt(2.0 * (kx**2) * cos2 + (ky**2) * sin2)
+    return np.sqrt(2.0 * ((kx**2) * cos2 + (ky**2) * sin2))
 
 
 def compute_E(kappa, c, L2D, z_i):
@@ -88,7 +88,7 @@ def generate_2D_lowfreq_approx(Nx, Ny, L1, L2, psi_degs, sigma2, L2D, z_i):
     if var_now > 1e-12:
         u_field *= np.sqrt(sigma2 / var_now)
 
-    return u_field
+    return np.linspace(0, L1, Nx), np.linspace(0, L2, Ny), u_field
 
 
 if __name__ == "__main__":

examples/08_mann_box_generation_IEC.py

@@ -27,9 +27,6 @@
 import numpy as np
 import torch
 
-from drdmannturb.fluctuation_generation import (
-    plot_velocity_components,  # utility function for plotting each velocity component in the field, not used in this example
-)
 from drdmannturb.fluctuation_generation import (
     GenerateFluctuationField,
     plot_velocity_magnitude,
@@ -97,9 +94,7 @@
     seed=seed,
 )
 
-fluctuation_field_mann = gen_mann.generate(
-    nBlocks, zref, uref, z0, windprofiletype, plexp
-)
+fluctuation_field_mann = gen_mann.generate(nBlocks, zref, uref, z0, windprofiletype, plexp)
 
 #######################################################################################
 # Adding the mean velocity profile
@@ -113,13 +108,11 @@
 
 spacing = tuple(grid_dimensions / (2.0**grid_levels + 1))
 
-fig_magnitude_mann = plot_velocity_magnitude(
-    spacing, fluctuation_field_mann, transparent=True
-)
+fig_magnitude_mann = plot_velocity_magnitude(spacing, fluctuation_field_mann, transparent=True)
 
-# this is a Plotly figure, which can be visualized with the ``.show()`` method in different contexts. While these utilities
-# may be useful for quick visualization, we recommend using Paraview to visualize higher resolution output. We will cover
-# saving to a portable VTK format further in this example.
+# this is a Plotly figure, which can be visualized with the ``.show()`` method in different contexts. While these
+# utilities may be useful for quick visualization, we recommend using Paraview to visualize higher resolution output.
+# We will cover saving to a portable VTK format further in this example.
 
 fig_magnitude_mann  # .show("browser"), or for specific browser, use .show("firefox")
 
@@ -141,10 +134,6 @@
 # -----------------------------------------
 # For higher resolution fluctuation fields, we suggest using Paraview. To transfer the generated data
 # from our package, we provide the ``.save_to_vtk()`` method.
-filename = str(
-    path / "./outputs/IEC_simple"
-    if path.name == "examples"
-    else path / "./outputs/IEC_simple"
-)
+filename = str(path / "./outputs/IEC_simple" if path.name == "examples" else path / "./outputs/IEC_simple")
 
 gen_mann.save_to_vtk(filename)