WIP solver restucture

docs/user_defined_wake_models.ipynb

@@ -0,0 +1,297 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "ba9ae6ce",
+   "metadata": {},
+   "source": [
+    "# User-defined Wake Models\n",
+    "\n",
+    "Beginning in v5, FLORIS supports user-defined wake models that can be passed directly into FLORIS using `fmodel.set_wake_model`. A user-defined wake model may be a dynamic or static class, but will usually be dynamic to allow model parameters to be set as attributes. It must conform to the `attrs` package for declaring attributes (in particular, wake model parameters). Additionally all user-defined operation models should inherit from the abstract parent class `BaseWakeModel`, available in FLORIS.\n",
+    "\n",
+    "All operation models must implement the following \"fundamental\" methods:\n",
+    "- `turbine_solve`: computes the flow solution at all turbine locations, part of the main FLORIS `run` procedure.\n",
+    "- `point_solve`: computes the flow solution at arbitrary, user-provided points in the flow, or for cut planes for visualization purposes.\n",
+    "\n",
+    "Wake models may then implement additional methods as needed.\n",
+    "\n",
+    "The following arguments are passed to either `turbine_solve` or `point_solve` at runtime:\n",
+    "\n",
+    "| Argument | Data type | Description |\n",
+    "|----------|-----------|----------|\n",
+    "| `farm` | `floris.core.Farm` | text |\n",
+    "| `flow_field` | `floris.core.FlowField` | The flow field object, which contains the flow solution and other flow-related quantities. |\n",
+    "| `grid` | `floris.core.TurbineGrid` or `floris.core.FlowFieldPlanarGrid` or `floris.core.PointsGrid` | The grid object corresponding to the type of solve being performed. For `turbine_solve`, this will be a `TurbineGrid`. For `point_solve`, this will be either a `FlowFieldPlanarGrid` or `PointsGrid`, depending on the type of points being solved for (visualization-type solves or individual point solves, respectively). |\n",
+    "\n",
+    "The `turbine_solve` and `point_solve` methods do not return any values, but instead update the `flow_field` argument in-place."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "aefe4f59",
+   "metadata": {},
+   "source": [
+    "### Static example\n",
+    "\n",
+    "We begin with a very simple example that will produce a \"straight\" wake behind each turbine, whose velocity deficit (as a fraction of the free stream velocity) is constant and user-definable. This is not a good wake model (and doesn't adhere to momentum conservation)! We're just using it as a basic example to demonstrate the functionality."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "d751aa3c",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from attrs import define, field\n",
+    "from floris.type_dec import floris_float_type, NDArrayFloat\n",
+    "from floris.core.wake_model import BaseWakeModel\n",
+    "from floris.flow_visualization import visualize_cut_plane\n",
+    "\n",
+    "from floris.core import (\n",
+    "    BaseModel,\n",
+    "    Farm,\n",
+    "    FlowField,\n",
+    "    FlowFieldPlanarGrid,\n",
+    "    PointsGrid,\n",
+    "    TurbineGrid,\n",
+    ")\n",
+    "\n",
+    "@define\n",
+    "class StraightWake(BaseWakeModel):\n",
+    "    \"\"\"\n",
+    "    A simple wake model that produces a straight wake behind each turbine.\n",
+    "    \"\"\"\n",
+    "\n",
+    "    # Using attrs, we can define model parameters as class attributes.\n",
+    "    velocity_deficit: float = field(default=0.2)\n",
+    "    wake_width: float = field(default=100.0)\n",
+    "\n",
+    "    # Define a method for determining whether a test point is within the wake of at least one\n",
+    "    # turbine\n",
+    "    def _is_in_wake(self, grid, turbine_i_x, turbine_i_y, turbine_i_z):\n",
+    "\n",
+    "        # Declare all True to start\n",
+    "        in_wake_i = np.full(grid.x_sorted.shape, True)\n",
+    "\n",
+    "        # Check if downstream of any turbine\n",
+    "        in_wake_i &= (grid.x_sorted > turbine_i_x.mean(axis=(2,3), keepdims=True))\n",
+    "\n",
+    "        # Check if within wake width of any turbine\n",
+    "        in_wake_i &= (\n",
+    "            np.abs(grid.y_sorted - turbine_i_y.mean(axis=(2,3), keepdims=True))\n",
+    "            < self.wake_width / 2\n",
+    "        )\n",
+    "        in_wake_i &= (\n",
+    "            np.abs(grid.z_sorted - turbine_i_z.mean(axis=(2,3), keepdims=True))\n",
+    "            < self.wake_width / 2\n",
+    "        )\n",
+    "\n",
+    "        # Return resulting boolean array\n",
+    "        return in_wake_i\n",
+    "\n",
+    "    # Define the main turbine_solve method for solving at turbine locations\n",
+    "    def turbine_solve(\n",
+    "        self,\n",
+    "        farm: Farm,\n",
+    "        flow_field: FlowField,\n",
+    "        grid: TurbineGrid,\n",
+    "    ) -> None:\n",
+    "\n",
+    "        # Initialize an array to keep track of whether each point is in the wake of any turbine\n",
+    "        in_wake = np.full(grid.x_sorted.shape, False)\n",
+    "\n",
+    "        for i in range(grid.n_turbines):\n",
+    "\n",
+    "            # Check if the points are in the wake of turbine i\n",
+    "            in_wake_i = self._is_in_wake(\n",
+    "                grid,\n",
+    "                grid.x_sorted[:, i:i+1, :, :],\n",
+    "                grid.y_sorted[:, i:i+1, :, :],\n",
+    "                grid.z_sorted[:, i:i+1, :, :]\n",
+    "            )\n",
+    "\n",
+    "            # Update the overall in_wake array to include the wake of turbine i\n",
+    "            in_wake |= in_wake_i\n",
+    "\n",
+    "        # Apply velocity deficits\n",
+    "        flow_field.u_sorted = flow_field.u_initial_sorted * (1 - self.velocity_deficit * in_wake)\n",
+    "\n",
+    "        return None\n",
+    "\n",
+    "    # Define the secondary point_solve method for solving at arbitrary points in the flow\n",
+    "    def point_solve(\n",
+    "        self,\n",
+    "        farm: Farm,\n",
+    "        flow_field: FlowField,\n",
+    "        grid: FlowFieldPlanarGrid | PointsGrid,\n",
+    "    ) -> None:\n",
+    "        # Use parent class method to access the turbine grid\n",
+    "        turbine_grid = self.generate_turbine_grid_objects(farm, flow_field)[2]\n",
+    "\n",
+    "        # Initialize an array to keep track of whether each point is in the wake of any turbine\n",
+    "        in_wake = np.full(grid.x_sorted.shape, False)\n",
+    "\n",
+    "        for i in range(turbine_grid.n_turbines):\n",
+    "\n",
+    "            # Check if the turbine is in the wake of any other turbine\n",
+    "            in_wake_i = self._is_in_wake(\n",
+    "                grid,\n",
+    "                turbine_grid.x_sorted[:, i:i+1, :, :],\n",
+    "                turbine_grid.y_sorted[:, i:i+1, :, :],\n",
+    "                turbine_grid.z_sorted[:, i:i+1, :, :]\n",
+    "            )\n",
+    "\n",
+    "            # Update the overall in_wake array to include the wake of turbine i\n",
+    "            in_wake |= in_wake_i\n",
+    "\n",
+    "        # Apply velocity deficits\n",
+    "        flow_field.u_sorted = flow_field.u_initial_sorted * (1 - self.velocity_deficit * in_wake)\n",
+    "\n",
+    "        return None"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "29ea6830",
+   "metadata": {},
+   "source": [
+    "Let's now use this straight wake model in FLORIS."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "b74802c2",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from floris import FlorisModel, TimeSeries\n",
+    "\n",
+    "fmodel = FlorisModel(\"defaults\")\n",
+    "time_series = TimeSeries(\n",
+    "    wind_directions=np.array([270.0, 270.0, 280.0]),\n",
+    "    wind_speeds=np.array([8.0, 10.0, 12.0]),\n",
+    "    turbulence_intensities=np.array([0.06, 0.06, 0.06]),\n",
+    ")\n",
+    "fmodel.set(\n",
+    "    layout_x = [0.0, 500.0],\n",
+    "    layout_y = [0.0, 0.0],\n",
+    "    wind_data=time_series,\n",
+    ")\n",
+    "fmodel.set_wake_model(StraightWake(velocity_deficit=0.2, wake_width=100.0))\n",
+    "\n",
+    "fmodel.run()\n",
+    "\n",
+    "print(\"Powers [W]:\\n\", fmodel.get_turbine_powers(), \"\\n\")\n",
+    "print(\"Thrust coefficients [-]:\\n\", fmodel.get_turbine_thrust_coefficients(), \"\\n\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "aa2b73b0",
+   "metadata": {},
+   "source": [
+    "## Visualization example\n",
+    "\n",
+    "Now, we will perform a flow visualization, which uses the `point_solve` method."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "93f1e99f",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fmodel = FlorisModel(\"defaults\")\n",
+    "fmodel.set(wind_speeds=[8.0], wind_directions=[280.0], turbulence_intensities=[0.06])\n",
+    "fmodel.set(\n",
+    "    layout_x = [0.0, 500.0],\n",
+    "    layout_y = [0.0, 0.0],\n",
+    ")\n",
+    "# TEMPORARY: reset the wake model (TODO: Work out how to avoid doing this)\n",
+    "fmodel.set_wake_model(StraightWake(velocity_deficit=0.2, wake_width=100.0))\n",
+    "horizontal_plane = fmodel.calculate_horizontal_plane(\n",
+    "    x_resolution=200,\n",
+    "    y_resolution=100,\n",
+    "    height=90.0,\n",
+    ")\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "visualize_cut_plane(\n",
+    "    horizontal_plane,\n",
+    "    ax=ax,\n",
+    "    label_contours=False,\n",
+    "    title=\"Horizontal Flow with Turbine Rotors and labels\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0f3409fc",
+   "metadata": {},
+   "source": [
+    "## Prepackaged wake models\n",
+    "\n",
+    "Naturally, prepackaged wake models can also be used in this way. Let's take a look at using the `Gauss` wake model from FLORIS, either as one of the preset defaults or by passing the class in directly."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "1dd9ed49",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from floris.core.wake_model import Gauss\n",
+    "\n",
+    "fmodel.set_wake_model(Gauss()) # Use Gauss defaults\n",
+    "horizontal_plane = fmodel.calculate_horizontal_plane(\n",
+    "    x_resolution=200,\n",
+    "    y_resolution=100,\n",
+    "    height=90.0,\n",
+    ")\n",
+    "\n",
+    "fig, ax = plt.subplots()\n",
+    "visualize_cut_plane(\n",
+    "    horizontal_plane,\n",
+    "    ax=ax,\n",
+    "    label_contours=False,\n",
+    "    title=\"Horizontal Flow with Turbine Rotors and labels\",\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "3ef42e63",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "floris",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.13.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

floris/core/__init__.py

@@ -22,7 +22,7 @@
 
 import floris.logging_manager
 
-from .base import BaseClass, BaseModel, State
+from .base import BaseClass, BaseLibrary, BaseModel, State
 from .turbine.turbine import (
     axial_induction,
     power,

floris/core/base.py

@@ -1,4 +1,5 @@
 
+import importlib
 from abc import abstractmethod
 from enum import Enum
 from typing import (
@@ -63,3 +64,33 @@ def prepare_function() -> dict:
     @abstractmethod
     def function() -> None:
         raise NotImplementedError("BaseModel.function")
+
+@define
+class BaseLibrary(BaseClass):
+    """
+    Base class that writes the name and module of the class into the attrs dictionary.
+    """
+    __classinfo__: dict = {"module": "", "name": ""}
+    def __attrs_post_init__(self) -> None:
+        #import ipdb; ipdb.set_trace()
+        self.__classinfo__ = {
+            "module": type(self).__module__,
+            "name": type(self).__name__
+        }
+
+    @staticmethod
+    def from_dict(data_dict):
+        """Recreate instance from dictionary with class information"""
+        if "__classinfo__" not in data_dict:
+            raise ValueError(
+                "Dictionary does not contain class information. ",
+                "Insure inheritance from BaseLibrary."
+            )
+        data_noinfo = data_dict.copy()
+        class_info = data_noinfo.pop("__classinfo__")
+
+        # Import the module and get the class
+        module = importlib.import_module(class_info["module"])
+        cls = getattr(module, class_info["name"])
+
+        return cls(**data_noinfo)