Merge pull request #60 from METHODS-Group/ndim-issue59

Non-dimensionalization with zref, ustar and pretty printing for issue 59

docs/source/uml_drd.rst

@@ -52,7 +52,6 @@ These interactive UML diagrams have an issue with rendering the correct arrow ty
             - sigma : float 
             - Uref : float 
             - zref : float 
-            - Iref : float 
             - domain : torch.Tensor
         }
         class LossParameters {

drdmannturb/fluctuation_generation/cmap_util.py

@@ -0,0 +1,49 @@
+"""
+Utilities and constants for colormaps used in plotting. 
+"""
+
+FIELD_COLORSCALE = [
+    [0.0, "#000059"],
+    [0.023809523809523808, "#0a1061"],
+    [0.047619047619047616, "#101e69"],
+    [0.07142857142857142, "#172f73"],
+    [0.09523809523809523, "#204380"],
+    [0.11904761904761904, "#29568a"],
+    [0.14285714285714285, "#336591"],
+    [0.16666666666666666, "#3d7499"],
+    [0.19047619047619047, "#4985a6"],
+    [0.21428571428571427, "#4985a6"],
+    [0.23809523809523808, "#5697b3"],
+    [0.26190476190476186, "#5ea1bb"],
+    [0.2857142857142857, "#63a7bf"],
+    [0.30952380952380953, "#77bcd1"],
+    [0.3333333333333333, "#92d1e0"],
+    [0.3571428571428571, "#a7dde8"],
+    [0.38095238095238093, "#c0eaf0"],
+    [0.40476190476190477, "#d2f4f7"],
+    [0.42857142857142855, "#d4e5e6"],
+    [0.45238095238095233, "#d5d8d6"],
+    [0.47619047619047616, "#d5d1ce"],
+    [0.5, "#d5c4be"],
+    [0.5238095238095237, "#d4b7af"],
+    [0.5476190476190476, "#d3b0a7"],
+    [0.5714285714285714, "#d2a99f"],
+    [0.5952380952380952, "#d2a69c"],
+    [0.6190476190476191, "#d1a499"],
+    [0.6428571428571428, "#d09f93"],
+    [0.6666666666666666, "#cb887a"],
+    [0.6904761904761905, "#c7796a"],
+    [0.7142857142857142, "#c06555"],
+    [0.738095238095238, "#b74d3e"],
+    [0.7619047619047619, "#ab2d24"],
+    [0.7857142857142857, "#a92b23"],
+    [0.8095238095238095, "#a62821"],
+    [0.8333333333333333, "#991818"],
+    [0.8571428571428571, "#8c1119"],
+    [0.8809523809523809, "#860f1d"],
+    [0.9047619047619047, "#800d20"],
+    [0.9285714285714285, "#730e2c"],
+    [0.9523809523809523, "#670e31"],
+    [0.9761904761904762, "#660e31"],
+    [1.0, "#591236"],
+]

drdmannturb/fluctuation_generation/fluctuation_field_generator.py

@@ -288,7 +288,55 @@ def _generate_block(self) -> np.ndarray:
 
         return wind
 
-    def generate(self, num_blocks: int) -> np.ndarray:
+    def _normalize_block(
+        self, curr_block: np.ndarray, roughness_height: float, friction_velocity: float
+    ) -> np.ndarray:
+        r"""Normalizes the generated field by the logarithmic profile
+
+        .. math ::
+
+            \left\langle U_1(z)\right\rangle=\frac{u_{\ast}}{\kappa} \ln \left(\frac{z}{z_0}+1\right)
+
+
+        where :math:`u_{\ast}` is the friction velocity and :math:`z_0` is the roughness height. More information on this normalization can be found in
+
+        J. JCSS, “Probabilistic model code,” Joint Committee on Structural Safety (2001).
+
+        Parameters
+        ----------
+        roughness_height : float
+            Roughness height :math:`z_0`.
+        friction_velocity : float
+            Ground friction velocity :math:`u_{\ast}`.
+
+        Returns
+        -------
+        np.ndarray
+            Fluctuation field normalized by the logarithmic profile.
+        """
+        if not np.any(curr_block):
+            raise ValueError(
+                "No fluctuation field has been generated, call the .generate() method first."
+            )
+
+        sd = np.sqrt(np.mean(curr_block**2))
+        curr_block /= sd
+
+        z_space = np.linspace(
+            0.0, self.grid_dimensions[2], 2 ** (self.grid_levels[2]) + 1
+        )
+        mean_profile_z = self.log_law(z_space, roughness_height, friction_velocity)
+
+        mean_profile = np.zeros_like(curr_block)
+        mean_profile[..., 0] = np.tile(
+            mean_profile_z.T, (mean_profile.shape[0], mean_profile.shape[1], 1)
+        )
+
+        return curr_block + mean_profile
+
+    def generate(
+        self, num_blocks: int, roughness_height: float, friction_velocity: float
+    ) -> np.ndarray:
         """Generates the full fluctuation field in blocks. The resulting field is stored as the ``total_fluctuation`` field of this object, allowing for all metadata of the object to be stored safely with the fluctuation field, and also reducing data duplication for post-processing; all operations can be performed on this public variable.
 
         .. warning::
@@ -312,8 +360,14 @@ def generate(self, num_blocks: int) -> np.ndarray:
             )
 
         for _ in range(num_blocks):
+            t_block = self._generate_block()
+
+            normed_block = self._normalize_block(
+                t_block, roughness_height, friction_velocity
+            )
+
             self.total_fluctuation = np.concatenate(
-                (self.total_fluctuation, self._generate_block()), axis=0
+                (self.total_fluctuation, normed_block), axis=0
             )
 
         return self.total_fluctuation
@@ -346,7 +400,7 @@ def normalize(
         """
         if not np.any(self.total_fluctuation):
             raise ValueError(
-                "No fluctuation field has been generated, call the .generate() method first!"
+                "No fluctuation field has been generated, call the .generate() method first."
             )
 
         sd = np.sqrt(np.mean(self.total_fluctuation**2))

drdmannturb/fluctuation_generation/wind_plot.py

@@ -6,6 +6,8 @@
 import plotly.graph_objects as go
 from plotly.subplots import make_subplots
 
+from .cmap_util import FIELD_COLORSCALE
+
 
 def create_grid(
     spacing: tuple[float, float, float], shape: tuple[int, int, int]
@@ -22,13 +24,13 @@ def create_grid(
     Returns
     -------
     np.ndarray
-       np.meshgrid object consisting of points at the provided spacing and with the specified counts in each dimension.
+       np.meshgrid object consisting of points at the provided spacing and with the specified counts in each dimension. This is 'ij' indexed (not Cartesian!).
     """
     x = np.array([spacing[0] * n for n in range(shape[0])])
     y = np.array([spacing[1] * n for n in range(shape[1])])
     z = np.array([spacing[2] * n for n in range(shape[2])])
 
-    return np.meshgrid(x, y, z)
+    return np.meshgrid(x, y, z, indexing="ij")
 
 
 def format_wind_field(
@@ -153,7 +155,7 @@ def plot_velocity_components(
         ),
     )
 
-    fig.update_coloraxes(colorscale="spectral_r")
+    fig.update_coloraxes(colorscale=FIELD_COLORSCALE)
 
     return fig
 
@@ -203,12 +205,21 @@ def plot_velocity_magnitude(
             z=Z.flatten(),
             value=wind_magnitude.flatten(),
             surface_count=surf_count,
-            opacity=0.75,
-            colorscale="Turbo",  # spectral_r
+            opacityscale=[
+                [0, 0.75],
+                [0.25, 0.5],
+                [0.35, 0.35],
+                [0.5, 0.1],
+                [0.75, 0.25],
+                [0.9, 0.35],
+                [1, 1],
+            ],
+            opacity=0.85,
+            colorscale=FIELD_COLORSCALE,
             colorbar={
                 "title": "|U(x)|",
             },
-            showscale=False,
+            showscale=True,
         ),
     )
 

drdmannturb/parameters.py

@@ -80,8 +80,6 @@ class PhysicalParameters:
         Reference velocity value at hub height (m/s)
     zref : float, optional
         Reference height value; should be measured at hub height (meters)
-    Iref : float, optional
-        Longitudinal turbulence scale parameter at hub height (meters)
     domain : torch.Tensor
         :math:`k_1` domain over which spectra data are defined.
     """
@@ -92,7 +90,8 @@ class PhysicalParameters:
 
     Uref: float = 10.0
     zref: float = 1.0
-    Iref: float = 0.14
+
+    ustar: float = 1.0
 
     domain: torch.Tensor = torch.logspace(-1, 2, 20)
 

drdmannturb/spectra_fitting/calibration.py

@@ -7,7 +7,7 @@
 from collections.abc import Iterable
 from functools import partial
 from pathlib import Path
-from typing import Optional, Union
+from typing import Dict, Optional, Union
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -101,6 +101,7 @@ def __init__(
 
         self.OPS = OnePointSpectra(
             type_eddy_lifetime=self.prob_params.eddy_lifetime,
+            physical_params=self.phys_params,
             type_power_spectra=self.prob_params.power_spectra,
             nn_parameters=self.nn_params,
             learn_nu=self.prob_params.learn_nu,
@@ -342,13 +343,15 @@ def calibrate(
         data: tuple[Iterable[float], torch.Tensor],
         tb_comment: str = "",
         optimizer_class: torch.optim.Optimizer = torch.optim.LBFGS,
-    ) -> np.ndarray:
+    ) -> Dict[str, float]:
         r"""Calibration method, which handles the main training loop and some
         data pre-processing. Currently the only supported optimizer is Torch's ``LBFGS``
         and the learning rate scheduler uses cosine annealing. Parameters for these
         components of the training process are set in ``LossParameters`` and ``ProblemParameters``
         during initialization.
 
+        See the ``.print_calibrated_params()`` method to receive a pretty-printed output of the calibrated physical parameters.
+
         Parameters
         ----------
         data : tuple[Iterable[float], torch.Tensor]
@@ -360,8 +363,8 @@ def calibrate(
 
         Returns
         -------
-        np.ndarray
-            Physical parameters for the problem, in normal space (not in log-space, as represented internally). In the order ``[length scale, time scale, spectrum amplitude]``.
+        Dict[str, float]
+            Physical parameters for the problem, in normal space (not in log-space, as represented internally). Presented as a dictionary in the order ``{L : length scale, Gamma : time scale, sigma : spectrum amplitude}``.
 
         Raises
         ------
@@ -382,9 +385,10 @@ def calibrate(
         self.k1_data_pts = torch.tensor(DataPoints, dtype=torch.float64)[:, 0].squeeze()
 
         self.LossAggregator = LossAggregator(
-            self.loss_params,
-            self.k1_data_pts,
-            self.logging_directory,
+            params=self.loss_params,
+            k1space=self.k1_data_pts,
+            zref=self.phys_params.zref,
+            tb_log_dir=self.logging_directory,
             tb_comment=tb_comment,
         )
 
@@ -479,12 +483,45 @@ def closure():
             print(f"Learned nu value: {self.OPS.tauNet.Ra.nu.item()}")
 
         # physical parameters are stored as natural logarithms internally
-        return np.exp(self.parameters)
+        self.calibrated_params = {
+            "L": np.exp(self.parameters[0]),
+            "Gamma": np.exp(self.parameters[1]),
+            "sigma": np.exp(self.parameters[2]),
+        }
+
+        return self.calibrated_params
 
     # ------------------------------------------------
     ### Post-treatment and Export
     # ------------------------------------------------
 
+    def print_calibrated_params(self):
+        """Prints out the optimal calibrated parameters ``L``, ``Gamma``, ``sigma``, which are stored in a fitted ``CalibrationProblem`` object under the ``calibrated_params`` dictionary.
+
+        Raises
+        ------
+        ValueError
+            Must call ``.calibrate()`` method to compute a fit to physical parameters first.
+        """
+        if not hasattr(self, "calibrated_params"):
+            raise ValueError(
+                "Must call .calibrate() method to compute a fit to physical parameters first."
+            )
+
+        OKGREEN = "\033[92m"
+        ENDC = "\033[0m"
+
+        print("=" * 40)
+
+        for k, v in self.calibrated_params.items():
+            print(
+                f"{OKGREEN}Optimal calibrated {k+' '*4 if k == 'L' else k} : {v:8.4f} {ENDC}"
+            )
+
+        print("=" * 40)
+
+        return
+
     def num_trainable_params(self) -> int:
         """Computes the number of trainable network parameters
             in the underlying model.
@@ -705,10 +742,14 @@ def plot(
                     "Must either provide data points or re-use what was used for model calibration, neither is currently specified."
                 )
 
-        kF_model_vals = model_vals if model_vals is not None else self.OPS(k1_data_pts)
+        kF_model_vals = (
+            model_vals
+            if model_vals is not None
+            else self.OPS(k1_data_pts) / self.phys_params.ustar**2.0
+        )
 
         kF_model_vals = kF_model_vals.cpu().detach()
-        kF_data_vals = kF_data_vals.cpu().detach()
+        kF_data_vals = kF_data_vals.cpu().detach() / self.phys_params.ustar**2
 
         if plot_tau:
             k_gd = torch.logspace(-3, 3, 50, dtype=torch.float64)

drdmannturb/spectra_fitting/data_generator.py

@@ -44,6 +44,7 @@ class OnePointSpectraDataGenerator:
     def __init__(
         self,
         zref: float,
+        ustar: float = 1.0,
         data_points: Optional[Iterable[Tuple[torch.tensor, float]]] = None,
         data_type: DataType = DataType.KAIMAL,
         k1_data_points: Optional[torch.Tensor] = None,
@@ -55,7 +56,9 @@ def __init__(
         Parameters
         ----------
         zref : float
-            Reference altitude value, by default 1.0
+            Reference altitude value
+        ustar : float
+            Friction velocity, by default 1.0.
         data_points : Iterable[Tuple[torch.tensor, float]], optional
             Observed spectra data points at each of the :math:`k_1` coordinates, paired with the associated reference height (typically kept at 1, but may depend on applications).
         data_type : DataType, optional
@@ -90,6 +93,7 @@ def __init__(
             self.spectra_values = spectra_values
 
         self.zref = zref
+        self.ustar = ustar
 
         if self.data_type == DataType.VK:
             self.eval = self.eval_VK
@@ -182,7 +186,7 @@ def generate_Initial_Parameters():
                 fit4 = fitOPS(self.k1, Data_temp[:, 3], 3)
                 DataValues[:, 0, 2] = -func3(self.k1, *fit4)
 
-                DataValues = torch.tensor(DataValues)
+                DataValues = torch.tensor(DataValues * self.ustar**2)
                 DataPoints = list(
                     zip(self.DataPoints, [self.zref] * len(self.DataPoints))
                 )
@@ -326,11 +330,12 @@ def eval_Kaimal(self, k1: float) -> torch.Tensor:
         """
         n = 1 / (2 * np.pi) * k1 * self.zref
         F = torch.zeros([3, 3])
-        F[0, 0] = 102 * n / (1 + 33 * n) ** (5 / 3)
-        F[1, 1] = 17 * n / (1 + 9.5 * n) ** (5 / 3)
-        F[2, 2] = 2.1 * n / (1 + 5.3 * n ** (5 / 3))
-        F[0, 2] = -12 * n / (1 + 9.6 * n) ** (7.0 / 3.0)
-        return F
+        F[0, 0] = 52.5 * n / (1 + 33 * n) ** (5 / 3)
+        F[1, 1] = 8.5 * n / (1 + 9.5 * n) ** (5 / 3)
+        F[2, 2] = 1.05 * n / (1 + 5.3 * n ** (5 / 3))
+        F[0, 2] = -7 * n / (1 + 9.6 * n) ** (12.0 / 5.0)
+
+        return F * self.ustar**2
 
     def plot(
         self,