Exclusion zone implementation

ard/api/interface.py

@@ -113,7 +113,7 @@ def set_up_ard_model(input_dict: Union[str, dict], root_data_path: str = None):
 def set_up_system_recursive(
     input_dict: dict,
     system_name: str = "top_level",
-    work_dir: str = "ard_prob_out",
+    work_dir: str = "case_files",
     parent_group=None,
     modeling_options: dict = None,
     analysis_options: dict = None,
@@ -132,7 +132,9 @@ def set_up_system_recursive(
     """
     # Initialize the top-level problem if no parent group is provided
     if parent_group is None:
-        prob = om.Problem(work_dir=work_dir)
+        prob = om.Problem(
+            work_dir=work_dir,
+        )
         parent_group = prob.model
         # parent_group.name = "ard_model"
     else:

ard/collection/optiwindnet_wrap.py

@@ -1,3 +1,5 @@
+from warnings import warn
+
 import networkx as nx
 import numpy as np
 
@@ -10,6 +12,7 @@
 
 
 def _own_L_from_inputs(inputs: dict, discrete_inputs: dict) -> nx.Graph:
+    # get the metadata and data for the OWN warm-starter from the inputs
     T = len(inputs["x_turbines"])
     R = len(inputs["x_substations"])
     name_case = "farm"
@@ -22,13 +25,63 @@ def _own_L_from_inputs(inputs: dict, discrete_inputs: dict) -> nx.Graph:
     VertexC[:T, 1] = inputs["y_turbines"]
     VertexC[-R:, 0] = inputs["x_substations"]
     VertexC[-R:, 1] = inputs["y_substations"]
+
+    # add perturbation to duplicate turbine/substation positions
+    VertexCTR = np.vstack([VertexC[:T, :], VertexC[-R:, :]])
+    perturbation_eps = 1.0e-6  # base magnitude of perturbation in m
+    perturbation_normal = np.array([-1.0, 1.0])  # set a fixed axis to perturb on
+    perturbation_normal = perturbation_normal / np.sqrt(
+        np.sum(perturbation_normal**2)
+    )  # normalize the perturbation
+    # go through the turbine/substation vertices and count how many times a
+    # given vertex has appeared before
+    repeat_accumulate = np.array(
+        [
+            int(np.sum(np.all(VertexCTR[:ivv, :] == vv, axis=1)))
+            for ivv, vv in enumerate(VertexCTR)
+        ]
+    )
+    if np.any(repeat_accumulate > 0):  # only if there are any repeats
+        warn_string = (
+            f"\nDetected {np.sum(repeat_accumulate > 0)} coincident "
+            f"turbines and/or substations in optiwindnet setup."
+        )  # start a warning string for the UserWarning
+        # TODO: make Ard warnings?
+
+        # create perturbation adjustements s.t. vertices w/ multiplicity > 2
+        # are adjusted to be fully unique!
+        adjustments = perturbation_eps * np.outer(
+            repeat_accumulate, perturbation_normal
+        )
+        # for each adjustments add to the warning string
+        for idx, dxy in enumerate(adjustments[:T, :]):
+            if np.sum(dxy != 0) == 0:
+                continue
+            warn_string += f"\n\tadjusting turbine #{idx} from {VertexCTR[idx, :]} to  {VertexCTR[idx, :] + dxy}"
+        for idx, dxy in enumerate((adjustments[-R:, :])[::-1, :]):
+            if np.sum(dxy != 0) == 0:
+                continue
+            warn_string += f"\n\tadjusting substation #{idx} from {VertexCTR[-(idx+1), :]} to {VertexCTR[-(idx+1), :] + dxy}"
+        # output the final warning
+        warn(warn_string)
+
+        # store the adjustments
+        VertexCTR += adjustments
+
+    # apply the adjustments
+    VertexC[:T, :] = VertexCTR[:T, :]
+    VertexC[-R:, :] = VertexCTR[-R:, :]
+
+    # put together the inputs for optiwindnet
     site = dict(
         T=T,
         R=R,
         name=name_case,
         handle=name_case,
         VertexC=VertexC,
     )
+
+    # handle the boundary if it exists
     if B > 0:
         VertexC[T:-R, 0] = discrete_inputs["x_border"]
         VertexC[T:-R, 1] = discrete_inputs["y_border"]

ard/eco/eagle_density.py

@@ -0,0 +1,84 @@
+import numpy as np
+
+import openmdao.api as om
+
+
+class EagleDensityFunction(om.ExplicitComponent):
+    """
+    _summary_
+
+    _extended_summary_
+
+    Options
+    -------
+    modeling_options : dict
+        a modeling options dictionary (inherited from
+        `templates.LanduseTemplate`)
+
+    Inputs
+    ------
+    x_turbines : np.ndarray
+        a 1-D numpy array that represents that x (i.e. Easting) coordinate of
+        the location of each of the turbines in the farm in meters
+    y_turbines : np.ndarray
+        a 1-D numpy array that represents that y (i.e. Northing) coordinate of
+        the location of each of the turbines in the farm in meters
+
+    Outputs
+    -------
+    area_tight : float
+        the area in square kilometers that the farm occupies based on the
+        circumscribing geometry with a specified (default zero) layback buffer
+        (inherited from `templates.LayoutTemplate`)
+    """
+
+    def initialize(self):
+        """Initialization of OM component."""
+        self.options.declare("modeling_options")
+
+    def setup(self):
+        """Setup of OM component."""
+
+        # load modeling options and turbine count
+        modeling_options = self.modeling_options = self.options["modeling_options"]
+        self.windIO = self.modeling_options["windIO_plant"]
+        self.N_turbines = modeling_options["layout"]["N_turbines"]
+
+        # self.eagle_density_function = lambda x, y: ???
+
+        # add the full layout inputs
+        self.add_input(
+            "x_turbines",
+            np.zeros((self.N_turbines,)),
+            units="m",
+            desc="turbine location in x-direction",
+        )
+        self.add_input(
+            "y_turbines",
+            np.zeros((self.N_turbines,)),
+            units="m",
+            desc="turbine location in y-direction",
+        )
+
+        # add outputs that are universal
+        self.add_output(
+            "eagle_normalized_density",
+            np.zeros((self.N_turbines,)),
+            units="unitless",
+            desc="normalized eagle presence density",
+        )
+
+
+    def compute(self, inputs, outputs):
+        """
+        Computation for the OM component.
+        """
+
+        x_turbines = inputs["x_turbines"]  # m
+        y_turbines = inputs["y_turbines"]  # m
+
+        raise NotImplementedError(
+            "@Eliot, you need to implement this!!!"
+        )
+
+        outputs["eagle_normalized_density"] = y  # on [0, 1]

ard/farm_aero/floris.py

@@ -6,6 +6,7 @@
 import floris
 import floris.turbine_library.turbine_utilities
 
+import ard.utils.logging as ard_logging
 import ard.farm_aero.templates as templates
 
 
@@ -190,6 +191,7 @@ def initialize(self):
         """Initialization-time FLORIS management."""
         self.options.declare("case_title")
 
+    @ard_logging.component_log_capture
     def setup(self):
         """Setup-time FLORIS management."""
 
@@ -340,21 +342,27 @@ def initialize(self):
         super().initialize()  # run super class script first!
         FLORISFarmComponent.initialize(self)  # FLORIS superclass
 
+    @ard_logging.component_log_capture
     def setup(self):
         super().setup()  # run super class script first!
         FLORISFarmComponent.setup(self)  # setup a FLORIS run
 
+    @ard_logging.component_log_capture
     def setup_partials(self):
         FLORISFarmComponent.setup_partials(self)
 
+    @ard_logging.component_log_capture
     def compute(self, inputs, outputs):
 
         # generate the list of conditions for evaluation
-        self.time_series = floris.TimeSeries(
-            wind_directions=np.degrees(np.array(self.wind_query.wind_directions)),
-            wind_speeds=np.array(self.wind_query.wind_speeds),
-            turbulence_intensities=np.array(self.wind_query.turbulence_intensities),
-        )
+        if type(self.wind_query) == floris.TimeSeries:
+            self.time_series = self.wind_query
+        else:
+            self.time_series = floris.TimeSeries(
+                wind_directions=np.degrees(np.array(self.wind_query.wind_directions)),
+                wind_speeds=np.array(self.wind_query.wind_speeds),
+                turbulence_intensities=np.array(self.wind_query.turbulence_intensities),
+            )
 
         # set up and run the floris model
         self.fmodel.set(
@@ -442,13 +450,16 @@ def initialize(self):
         super().initialize()  # run super class script first!
         FLORISFarmComponent.initialize(self)  # add on FLORIS superclass
 
+    @ard_logging.component_log_capture
     def setup(self):
         super().setup()  # run super class script first!
         FLORISFarmComponent.setup(self)  # setup a FLORIS run
 
+    @ard_logging.component_log_capture
     def setup_partials(self):
         super().setup_partials()
 
+    @ard_logging.component_log_capture
     def compute(self, inputs, outputs):
 
         # set up and run the floris model
@@ -477,5 +488,6 @@ def compute(self, inputs, outputs):
         outputs["power_turbines"] = FLORISFarmComponent.get_power_turbines(self)
         outputs["thrust_turbines"] = FLORISFarmComponent.get_thrust_turbines(self)
 
+    @ard_logging.component_log_capture
     def setup_partials(self):
         FLORISFarmComponent.setup_partials(self)

ard/farm_aero/templates.py

@@ -10,6 +10,7 @@
 def create_windresource_from_windIO(
     windIOdict: dict,
     resource_type: str = None,  # ["probability", "timeseries", "weibull_sector"]
+    heterogeneous_map: floris.HeterogeneousMap | None = None,
 ):
     """
     takes a windIO plant specification and creates an appropriate wind resource
@@ -120,6 +121,7 @@ def create_windresource_from_windIO(
             wind_speeds=wind_speeds,
             freq_table=probabilities,
             ti_table=turbulence_intensities,
+            heterogeneous_map=heterogeneous_map,
         )
         # stash some metadata for the wind resource
         wind_resource_representation.reference_height = (
@@ -163,6 +165,7 @@ def create_windresource_from_windIO(
             wind_directions=wind_directions,
             wind_speeds=wind_speeds,
             turbulence_intensities=turbulence_intensities,
+            heterogeneous_map=heterogeneous_map,
         )
         # stash some metadata for the wind resource
         wind_resource_representation.reference_height = (
@@ -300,9 +303,25 @@ def setup(self):
         super().setup()
 
         # unpack wind query object
+        heterogeneous_map_spec = self.modeling_options.get("heterogeneous_map")
+        if heterogeneous_map_spec:
+            print("ACTIVATING HETEROGENEOUS MAP SPEC")
         self.wind_query = create_windresource_from_windIO(
             self.windIO,
             "timeseries",
+            heterogeneous_map=(
+                floris.HeterogeneousMap(
+                    x=heterogeneous_map_spec["x"],
+                    y=heterogeneous_map_spec["y"],
+                    z=heterogeneous_map_spec.get("z"),
+                    speed_multipliers=heterogeneous_map_spec["speed_multipliers"],
+                    wind_directions=heterogeneous_map_spec["wind_directions"],
+                    wind_speeds=heterogeneous_map_spec["wind_speeds"],
+                    interp_method=heterogeneous_map_spec.get("interp_method"),
+                )
+                if heterogeneous_map_spec
+                else None
+            ),
         )
         self.directions_wind = self.wind_query.wind_directions.tolist()
         self.speeds_wind = self.wind_query.wind_speeds.tolist()
@@ -428,9 +447,25 @@ def setup(self):
         super().setup()
         data_path = str(self.options["data_path"])
 
+        heterogeneous_map_spec = self.modeling_options.get("heterogeneous_map")
+        if heterogeneous_map_spec:
+            print("ACTIVATING HETEROGENEOUS MAP SPEC")
         self.wind_query = create_windresource_from_windIO(
             self.windIO,
             "probability",
+            heterogeneous_map=(
+                floris.HeterogeneousMap(
+                    x=heterogeneous_map_spec["x"],
+                    y=heterogeneous_map_spec["y"],
+                    z=heterogeneous_map_spec.get("z"),
+                    speed_multipliers=heterogeneous_map_spec["speed_multipliers"],
+                    wind_directions=heterogeneous_map_spec["wind_directions"],
+                    wind_speeds=heterogeneous_map_spec["wind_speeds"],
+                    interp_method=heterogeneous_map_spec.get("interp_method"),
+                )
+                if heterogeneous_map_spec
+                else None
+            ),
         )
 
         if data_path is None: