py-wake - npm Package Compare versions

Comparing version 2.6.13 to 2.6.14

py_wake/site/streamline_distance.py

+import warnings
+
+import matplotlib.pyplot as plt
+from numpy import newaxis as na
+
+from py_wake import np
+from py_wake.deficit_models.noj import NOJ
+from py_wake.examples.data.hornsrev1 import V80
+from py_wake.examples.data.ParqueFicticio._parque_ficticio import ParqueFicticioSite
+from py_wake.flow_map import XYGrid
+from py_wake.site.distance import StraightDistance
+from py_wake.utils.model_utils import DeprecatedModel
+from py_wake.utils.streamline import VectorField3D
+
+
+class StreamlineDistance(StraightDistance):
+    """Just-In-Time Streamline Distance
+    Calculates downwind crosswind and vertical distance along streamlines.
+    Streamlines calculated in each call
+    """
+
+    def __init__(self, vectorField, step_size=20):
+        """Parameters
+        ----------
+        vectorField : VectorField3d
+        step_size : int for float
+            Size of linear streamline steps
+        """
+        StraightDistance.__init__(self, wind_direction='wd')
+        self.vectorField = vectorField
+        self.step_size = step_size
+
+    def __call__(self, src_x_ilk, src_y_ilk, src_h_ilk, wd_l=None, WD_ilk=None, time=None, dst_xyh_jlk=None):
+        (src_x_ilk, src_y_ilk, src_h_ilk), (dst_x_jlk, dst_y_jlk, dst_h_jlk) = self.get_pos(
+            src_x_ilk, src_y_ilk, src_h_ilk, wd_l, WD_ilk, dst_xyh_jlk)
+        assert src_x_ilk.shape[2] == 1, 'StreamlineDistance does not support flowcase dependent positions'
+
+        start_points_m = np.moveaxis([v[:, :, 0].flatten() for v in [src_x_ilk, src_y_ilk, src_h_ilk]], 0, -1)
+
+        dw_ijlk, hcw_ijlk, dh_ijlk = StraightDistance.__call__(self, src_x_ilk, src_y_ilk, src_h_ilk, wd_l=wd_l,
+                                                               dst_xyh_jlk=dst_xyh_jlk)
+        src_z_ilk = self.site.elevation(src_x_ilk, src_y_ilk)
+        dst_z_jlk = self.site.elevation(dst_x_jlk, dst_y_jlk)
+        dz_ijlk = dst_z_jlk[na, :] - src_z_ilk[:, na]
+
+        # +0 ~ autograd safe copy (broadcast_to returns readonly array)
+        dh_ijlk = np.broadcast_to(dh_ijlk, dw_ijlk.shape) + 0.
+        dz_ijlk = np.broadcast_to(dz_ijlk, dw_ijlk.shape) + 0.
+        I, J, L, K = dw_ijlk.shape
+        dw_mj, hcw_mj, dh_mj, dz_mj = [np.moveaxis(v, 1, 2).reshape(I * L, J)
+                                       for v in [dw_ijlk, hcw_ijlk, dh_ijlk, dz_ijlk]]
+
+        wd_m = np.tile(wd_l, I)
+
+        stream_lines = self.vectorField.stream_lines(wd_m, time=time, start_points=start_points_m, dw_stop=dw_mj.max(1),
+                                                     step_size=self.step_size)
+
+        dxyz = np.diff(np.concatenate([stream_lines[:, :1], stream_lines], 1), 1, -2)
+        length_is = np.cumsum(np.sqrt(np.sum(dxyz**2, -1)), -1)
+        dist_xyz = stream_lines - start_points_m[:, na]
+        t = np.deg2rad(270 - wd_m)[:, na]
+        dw_is = dist_xyz[:, :, 0] * np.cos(t) + dist_xyz[:, :, 1] * np.sin(t)
+        hcw_is = dist_xyz[:, :, 0] * np.sin(t) - dist_xyz[:, :, 1] * np.cos(t)
+
+        for m, (dw_j, dw_s, hcw_s, dh_s, length_s) in enumerate(
+                zip(dw_mj, dw_is, hcw_is, dist_xyz[:, :, 2], length_is)):
+            dw = dw_j > 0
+            hcw_mj[m, dw] += np.interp(dw_j[dw], dw_s, hcw_s)
+            dh_mj[m, dw] -= np.interp(dw_j[dw], dw_s, dh_s)
+            dw_mj[m, dw] = np.interp(dw_j[dw], dw_s, length_s)
+
+        # streamline dh contains absolute height different, but pywake needs differences relative to ground, so
+        # we need to subtract elevation differences, dz
+        dh_mj += dz_mj
+
+        return [np.moveaxis(v.reshape((I, L, J, 1)), 2, 1) for v in [dw_mj, hcw_mj, dh_mj]]
+
+
+def main():
+    if __name__ == '__main__':
+
+        wt = V80()
+        vf3d = VectorField3D.from_WaspGridSite(ParqueFicticioSite())
+        site = ParqueFicticioSite(distance=StreamlineDistance(vf3d))
+
+        x, y = site.initial_position[:].T
+        wfm = NOJ(site, wt)
+        wd = 330
+        sim_res = wfm(x, y, wd=[wd], ws=10)
+        fm = sim_res.flow_map(XYGrid(x=np.linspace(site.ds.x[0].item(), site.ds.x[-1].item(), 500),
+                                     y=np.linspace(site.ds.y[0].item(), site.ds.y[-1].item(), 500)))
+        stream_lines = vf3d.stream_lines(wd=np.full(x.shape, wd), start_points=np.array([x, y, np.full(x.shape, 70)]).T,
+                                         dw_stop=y - 6504700)
+        fm.plot_wake_map()
+        for sl in stream_lines:
+            plt.plot(sl[:, 0], sl[:, 1])
+
+        plt.show()
+
+
+main()

py_wake/tests/test..._farm_models/test_external_windfarms.py

+import os
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+import xarray as xr
+from numpy import newaxis as na
+
+from py_wake.deficit_models import BastankhahGaussianDeficit, SelfSimilarityDeficit2020
+from py_wake.examples.data.hornsrev1 import V80, Hornsrev1Site
+from py_wake.flow_map import Points, XYGrid
+from py_wake.site.streamline_distance import StreamlineDistance
+from py_wake.superposition_models import LinearSum
+from py_wake.tests import npt
+from py_wake.utils import layouts
+from py_wake.utils.plotting import setup_plot
+from py_wake.utils.profiling import profileit, timeit
+from py_wake.utils.streamline import VectorField3D
+from py_wake.wind_farm_models import All2AllIterative, PropagateDownwind
+from py_wake.wind_farm_models.external_wind_farm_models import (
+    ExternalWFMWindFarm,
+    ExternalWindFarm,
+    ExternalXRAbsWindFarm,
+    ExternalXRRelWindFarm,
+)
+from py_wake.wind_turbines import WindTurbines
+
+
+def get_wfm(externalWindFarms=[], wfm_cls=PropagateDownwind, site=Hornsrev1Site(), blockage=False):
+    windTurbines = WindTurbines.from_WindTurbine_lst([V80()] * 6)
+    windTurbines._names = ["Current WF"] + [f"WF{i + 1}" for i in np.arange(5)]
+    kwargs = dict(site=site, windTurbines=windTurbines,
+                  wake_deficitModel=BastankhahGaussianDeficit(use_effective_ws=True),
+                  superpositionModel=LinearSum(),
+                  externalWindFarms=externalWindFarms)
+    if blockage:
+        wfm_cls = All2AllIterative
+        kwargs['blockage_deficitModel'] = SelfSimilarityDeficit2020()
+    return wfm_cls(**kwargs)
+
+
+def setup_ext_farms(cls, neighbour_x_y_angle, wfm=get_wfm(), include_wd_range=np.arange(-45, 45)):
+    ext_farms = []
+    for i, (x, y, angle) in enumerate(neighbour_x_y_angle, 1):
+        wd_lst = (include_wd_range + angle) % 360  # relevant wind directions
+        name = f'WF{i}'
+        if cls is ExternalWFMWindFarm:
+            ext_wf = ExternalWFMWindFarm(name, wfm, x, y, include_wd=wd_lst)
+        elif cls is ExternalXRRelWindFarm:
+            # Coarse grid in relative downwind, crosswind and vertical direction
+            grid_xyh = (np.linspace(1200, 2500, 12),
+                        np.linspace(-1000, 1000, 20),
+                        np.array([0]))
+            ext_wf = ExternalXRRelWindFarm.generate(name, grid_xyh, wfm, wt_x=x, wt_y=y, wd=wd_lst)
+        elif cls is ExternalXRAbsWindFarm:
+            # coarse grid in East, North and vertical direction covering the current wind farm
+            e = 500
+            grid_xyh = (np.linspace(- e, e, 20),
+                        np.linspace(- e, e, 20),
+                        wfm.windTurbines.hub_height())
+            ext_wf = ExternalXRAbsWindFarm.generate(name, grid_xyh, wfm, wt_x=x, wt_y=y, wd=wd_lst)
+        ext_farms.append(ext_wf)
+    return ext_farms
+
+
+def test_aep():
+    # setup current, neighbour and all positions
+    wf_x, wf_y = layouts.circular([1, 5], 400)
+
+    wts = WindTurbines.from_WindTurbine_lst([V80()] * 6)
+    wts._names = ["Current WF"] + [f"WF{i + 1}" for i in np.arange(5)]
+
+    No_neighbours = 3
+    neighbour_x_y_angle = [(wf_x + 2000 * np.cos(d), wf_y + 2000 * np.sin(d), (90 - np.rad2deg(d)) % 360)
+                           for d in np.pi + np.linspace(0, np.pi / 2, No_neighbours)]
+    neighbour_x, neighbour_y, _ = zip(*neighbour_x_y_angle)
+
+    all_x, all_y = np.r_[wf_x, np.array(neighbour_x).flatten()], np.r_[wf_y, np.array(neighbour_y).flatten()]
+
+    for wfm_cls in [PropagateDownwind, All2AllIterative]:
+        def run():
+            return get_wfm(wfm_cls=wfm_cls)(all_x, all_y, type=0).aep().isel(wt=np.arange(len(wf_x)))
+        aep_ref, t = timeit(run)()
+        # print(wfm_cls.__name__, np.mean(t), aep_ref.sum().item())
+
+        ext_cls_lst = [ExternalWFMWindFarm, ExternalXRAbsWindFarm, ExternalXRRelWindFarm]
+        for cls in ext_cls_lst:
+            def run():
+                ext_farm = setup_ext_farms(cls, neighbour_x_y_angle)
+                wfm = get_wfm(ext_farm, wfm_cls=wfm_cls)
+                if 0:
+                    wfm(wf_x, wf_y, type=0, wd=254, ws=10).flow_map().plot_wake_map()
+                    plt.title(f'{wfm_cls.__name__}, {cls.__name__}')
+                    plt.show()
+                return wfm(wf_x, wf_y, type=0).aep().isel(wt=np.arange(len(wf_x)))
+            aep, t = timeit(run, )()
+
+            err_msg = f'{wfm_cls.__name__}, {cls.__name__}'
+            atol = {ExternalXRAbsWindFarm: 0.0006, ExternalXRRelWindFarm: 0.001}.get(cls, 1e-6)
+            # print(err_msg, np.mean(t), aep.sum().item())
+            npt.assert_allclose(aep.sum().item(), aep_ref.sum().item(), atol=atol, err_msg=err_msg)
+            npt.assert_allclose(aep.values, aep_ref.values, atol=atol, err_msg=err_msg)
+
+
+def test_functionality():
+    wf_x, wf_y = layouts.circular([1, 5], 400)
+
+    wts = WindTurbines.from_WindTurbine_lst([V80()] * 6)
+    wts._names = ["Current WF"] + [f"WF{i + 1}" for i in np.arange(5)]
+
+    No_neighbours = 3
+    neighbour_x_y_angle = [(wf_x + 2000 * np.cos(d), wf_y + 2000 * np.sin(d), (90 - np.rad2deg(d)) % 360)
+                           for d in np.pi + np.linspace(0, np.pi / 2, No_neighbours)]
+    neighbour_x, neighbour_y, _ = zip(*neighbour_x_y_angle)
+
+    all_x, all_y = np.r_[wf_x, np.array(neighbour_x).flatten()], np.r_[wf_y, np.array(neighbour_y).flatten()]
+    ext_farm = setup_ext_farms(ExternalWFMWindFarm, neighbour_x_y_angle)
+    wfm = get_wfm(ext_farm)
+    wfm(wf_x, wf_y, wd=234, ws=10).flow_map().plot_wake_map()
+    if 0:
+        plt.show()
+    plt.close('all')
+
+    npt.assert_array_equal([wd for wd in np.arange(360) if ext_farm[0].include_wd_func(wd)], np.arange(225, 315))
+    ext_farm[0].set_include_wd(np.arange(225, 300))
+    npt.assert_array_equal([wd for wd in np.arange(360) if ext_farm[0].include_wd_func(wd)], np.arange(225, 300))
+    ext_farm[0].set_include_wd(lambda wd: 235 <= wd <= 245)
+    npt.assert_array_equal([wd for wd in np.arange(360) if ext_farm[0].include_wd_func(wd)], np.arange(235, 246))
+
+    npt.assert_array_equal(ext_farm[0].get_relevant_wd((wf_x, wf_y, wts.hub_height())), np.arange(236, 305))
+    npt.assert_array_equal(ext_farm[0].get_relevant_wd((wf_x, wf_y, wts.hub_height()), ws=18), np.arange(237, 304))
+    npt.assert_array_equal(ext_farm[0].get_relevant_wd((wf_x, wf_y, wts.hub_height()), tol=1e-3), np.arange(240, 301))
+
+    ext_farm = setup_ext_farms(ExternalXRAbsWindFarm, neighbour_x_y_angle[:1])
+    assert ext_farm[0].to_xarray().deficit.shape == (20, 20, 1, 90, 23)
+
+
+def test_cluster_interaction():
+    wf_x, wf_y = layouts.circular([1, 5, 12, 18], 1800)
+    No_neighbours = 2
+    neighbour_x_y_angle = [(wf_x - 6000 * i, wf_y, 270) for i in range(1, No_neighbours + 1)]
+    neighbour_x, neighbour_y, _ = zip(*neighbour_x_y_angle)
+    all_x, all_y = np.r_[wf_x, np.array(neighbour_x).flatten()], np.r_[wf_y, np.array(neighbour_y).flatten()]
+
+    types = [v for i in range(No_neighbours + 1) for v in [i] * len(wf_x)]
+    wd = 270
+    wfm_ref = get_wfm()
+    sim_res_ref = wfm_ref(all_x, all_y, type=types, wd=wd, ws=10)
+
+    ext_farms = setup_ext_farms(ExternalWFMWindFarm, neighbour_x_y_angle)
+    wfm = get_wfm(ext_farms)
+    sim_res_ext = wfm(wf_x, wf_y, type=0, wd=wd, ws=10)
+    if 0:
+        grid = XYGrid(x=np.linspace(-15000, 3000, 150), y=np.linspace(-2500, 2500, 100))
+        fm_ref = sim_res_ref.flow_map(grid)
+        fm_ref.plot_wake_map(levels=np.linspace(4, 10, 50))
+        setup_plot(grid=False, figsize=(12, 3), axis='scaled')
+        plt.show()
+    npt.assert_allclose(sim_res_ext.Power[:len(wf_x)].sum().item(),
+                        sim_res_ref.Power[:len(wf_x)].sum().item(), atol=26000)
+
+
+def test_cluster_blockage():
+    wf_x, wf_y = layouts.circular([1, 5, 12, 18], 1800)
+
+    No_neighbours = 1
+    neighbour_x_y_angle = [(wf_x - 6000 * i, wf_y, 270) for i in range(1, No_neighbours + 1)]
+    neighbour_x, neighbour_y, _ = zip(*neighbour_x_y_angle)
+    all_x, all_y = np.r_[wf_x, np.array(neighbour_x).flatten()], np.r_[wf_y, np.array(neighbour_y).flatten()]
+    types = [v for i in range(No_neighbours + 1) for v in [i] * len(wf_x)]
+
+    wfm_ref = get_wfm(blockage=True)
+    sim_res_ref = wfm_ref(all_x, all_y, type=types, wd=270, ws=10)
+
+    ext_farms = setup_ext_farms(ExternalWFMWindFarm, neighbour_x_y_angle, wfm=wfm_ref)
+    wfm = get_wfm(ext_farms, blockage=True)
+    sim_res_ext = wfm(wf_x, wf_y, type=0, wd=270, ws=10)
+
+    # all WT without blockage
+    npt.assert_allclose(sim_res_ref.Power[:len(wf_x)].sum().item(),
+                        get_wfm()(all_x, all_y, wd=270, ws=10).Power[:len(wf_x)].sum().item(), atol=9000)
+    # Current + external farms with blockage
+    npt.assert_allclose(sim_res_ref.Power[:len(wf_x)].sum().item(),
+                        sim_res_ext.Power[:len(wf_x)].sum().item(), atol=400)
+
+
+def test_streamlines():
+    wf_x, wf_y = layouts.circular([1, 5], 400)
+    H = 70
+
+    class MyVectorField(VectorField3D):
+        def __init__(self):
+            pass
+
+        def __call__(self, wd, time, x, y, h):
+            turning = (x + 1000) / 50
+            theta = np.deg2rad(270 - wd + turning)
+            return np.array([np.cos(theta), np.sin(theta), theta * 0]).T
+
+    vf3d = MyVectorField()
+    site = Hornsrev1Site()
+    site_with_streamlines = Hornsrev1Site()
+
+    site_with_streamlines.distance = StreamlineDistance(vf3d)
+
+    wfm_ref = get_wfm(site=site_with_streamlines)
+    No_neighbours = 1
+    neighbour_x_y_angle = (wf_x - 2000, wf_y, 270)
+    neighbour_x, neighbour_y, _ = neighbour_x_y_angle
+    all_x, all_y = np.r_[wf_x, np.array(neighbour_x).flatten()], np.r_[wf_y, np.array(neighbour_y).flatten()]
+    types = [v for i in range(No_neighbours + 1) for v in [i] * len(wf_x)]
+
+    grid = XYGrid(x=np.linspace(-4000, 2000, 150), y=np.linspace(-1500, 1500, 200))
+
+    sim_res_ref = wfm_ref(all_x, all_y, type=types, wd=270, ws=10)
+    df = pd.DataFrame({'Power current WF [MW]': []})
+    P_ref = sim_res_ref.Power[:len(wf_x)].sum().item()
+
+    if 0:
+        fm_ref = sim_res_ref.flow_map(grid)
+        fm_ref.plot_wake_map()
+        stream_lines = vf3d.stream_lines(wd=np.full_like(neighbour_x, 270), start_points=np.array([neighbour_x, neighbour_y, np.full_like(neighbour_x, 70)]).T,
+                                         dw_stop=np.full_like(neighbour_x, 2000))
+        for sl in stream_lines:
+            plt.plot(sl[:, 0], sl[:, 1], 'k', lw=1, alpha=0.1)
+
+        setup_plot(axis='scaled', xlabel='', ylabel='', xlim=[-4000, 2000], grid=0, figsize=(12, 4))
+        plt.show()
+
+    for cls in [ExternalWFMWindFarm, ExternalXRAbsWindFarm, ExternalXRRelWindFarm]:
+        fig, axes = plt.subplots(2, 1, figsize=(8, 6))
+
+        for i, (ax, site_setup) in enumerate(zip(axes, [site_with_streamlines, site]), 1):
+            # setup external wind farm model
+            # - option 1: with streamlines
+            # - option 2: without streamlines
+            ext_farms = setup_ext_farms(cls, [neighbour_x_y_angle], wfm=get_wfm(site=site_setup),
+                                        include_wd_range=np.array([0]))
+
+            # Setup simulation using site with stream lines
+            wfm = get_wfm(ext_farms, site=site_with_streamlines)
+            sim_res_ext = wfm(wf_x, wf_y, type=0, wd=270, ws=10)
+            s = f'{cls.__name__}, option {i}'
+            if i == 2:
+                rtol = 0.05
+            else:
+                rtol = {ExternalWFMWindFarm: 1e-7,
+                        ExternalXRAbsWindFarm: 0.0005,
+                        ExternalXRRelWindFarm: 0.009}[cls]
+            npt.assert_allclose(sim_res_ext.Power[:len(wf_x)].sum().item(), P_ref, rtol=rtol, err_msg=s)
+
+            if 0:
+                fm = sim_res_ext.flow_map(grid)
+
+                fm.plot_wake_map(ax=ax)
+
+                if site_setup == site_with_streamlines:
+                    start_points, alpha = np.array(
+                        [ext_farms[0].wt_x, ext_farms[0].wt_y, ext_farms[0].wt_x * 0 + H]).T, 0.2
+                else:
+                    start_points, alpha = np.array([[ext_farms[0].wf_x, ext_farms[0].wf_y, 70]]), .5
+                stream_lines = vf3d.stream_lines(wd=[270], start_points=start_points,
+                                                 dw_stop=np.full(start_points[:, 0].shape, 2000))
+                for sl in stream_lines:
+                    ax.plot(sl[:, 0], sl[:, 1], alpha=alpha)
+
+                setup_plot(ax=ax, axis='scaled', xlabel='', ylabel='', title=f'{cls.__name__}, option {i}', grid=0)
+
+    # plt.show()
+    plt.close('all')

py_wake/wind_farm_models/external_wind_farm_models.py

+from abc import ABC, abstractmethod
+
+import matplotlib.pyplot as plt
+import numpy as np
+import xarray as xr
+from numpy import newaxis as na
+from scipy.interpolate import RegularGridInterpolator as RGI
+from tqdm import tqdm
+
+from py_wake.flow_map import Points
+from py_wake.site.distance import StraightDistance
+from py_wake.utils.gradients import item_assign
+from py_wake.utils.model_utils import Model
+from py_wake.utils.xarray_utils import sel_interp_all
+
+
+class ExternalWindFarm(Model, ABC):
+
+    def __init__(self, name, windTurbines, wt_x, wt_y, wt_h=None):
+        self.name = name
+        self.wt_x = wt_x
+        self.wt_y = wt_y
+        self.wf_x = np.mean([np.min(wt_x), np.max(wt_x)])
+        self.wf_y = np.mean([np.min(wt_y), np.max(wt_y)])
+        if wt_h is None:
+            self.wt_h = windTurbines.hub_height(windTurbines.types())
+        self.wf_h = np.mean(self.wt_h)
+        self.windTurbines = windTurbines
+
+    def set_include_wd(self, include_wd):
+        if callable(include_wd):
+            self.include_wd_func = include_wd
+        else:
+            assert isinstance(include_wd, (int, float, list, tuple, np.ndarray))
+
+            def include_wd_func(wd):
+                return np.round(wd) in set(np.asarray(include_wd).tolist())
+            self.include_wd_func = include_wd_func
+
+    def rel2abs(self, dw, hcw, dh, WD):
+        theta = np.deg2rad(270 - WD)
+        co, si = np.cos(theta), np.sin(theta)
+        x = co * dw - hcw * si + self.wf_x
+        y = si * dw + hcw * co + self.wf_y
+        h = dh + self.wf_h
+        return x, y, h
+
+    def plot(self, ax=None):
+        ax = ax or plt.gca()
+        ax.plot(self.wt_x, self.wt_y, '1', label=self.name)
+
+    @abstractmethod
+    def __call__(self):
+        ""
+
+
+class ExternalWFMWindFarm(ExternalWindFarm):
+    def __init__(self, name, windFarmModel, wt_x, wt_y, wt_h=None, type=0, include_wd=np.arange(360), **kwargs):
+        self.wfm = windFarmModel
+        self.wfm_kwargs = {**kwargs, 'x': wt_x, 'y': wt_y, 'h': wt_h, 'type': type}
+        self.set_include_wd(include_wd)
+        ExternalWindFarm.__init__(self, name, windFarmModel.windTurbines, wt_x, wt_y, wt_h)
+
+    def get_relevant_wd(self, target_xyh, ws=10, tol=1e-6):
+        x, y, h = target_xyh
+        h = np.zeros_like(x) + h
+        sim_res = self.get_sim_res(wd=np.arange(360), ws=[ws])
+        fm = sim_res.flow_map(Points(x, y, h), wd=sim_res.wd)
+        return fm.wd[((fm.WS - fm.WS_eff).max('i') > tol)].values
+
+    def get_sim_res(self, wd, ws):
+        return self.wfm(**{**self.wfm_kwargs, 'ws': ws, 'wd': wd})
+
+    def __call__(self, i, l, deficit_jlk, WS_jlk, WS_eff_ilk,
+                 WD_ilk, dst_xyh_jlk, IJLK, dw_ijlk, hcw_ijlk, dh_ijlk, **_):
+        I, J, L, K = IJLK
+        WD_l = np.round(WD_ilk[np.minimum(i, len(WD_ilk) - 1), :, 0])
+        m_lst = [m for m, wd in enumerate(WD_l) if l[m] and self.include_wd_func(wd)]
+        if deficit_jlk is None:
+            deficit_jlk = np.zeros((J, L, K))
+        if m_lst:
+            M = len(m_lst)
+            ws_lst = np.sort(np.unique(np.r_[WS_eff_ilk.min(), np.round(WS_eff_ilk.flatten()), WS_eff_ilk.max()]))
+            wd_lst = WD_l[m_lst]
+            sim_res = self.wfm(**{**self.wfm_kwargs, 'ws': ws_lst, 'wd': np.sort(np.unique(wd_lst))})
+            from py_wake.site.streamline_distance import StreamlineDistance
+            if isinstance(self.wfm.site.distance, StreamlineDistance):
+                # Simulation and flow map of external WF is computed with stream line distances
+                # the WF wake must therefore be calculated for fixed abs positions
+                dst_xyh_jlk = dst_xyh_jlk
+            else:
+                # Simulation and flow map of external WF is computed with straight distance.
+                # The WF wakes must therefore be calcualted using the relative distances (including streamlines)
+                _dst_xyh_jlk = self.rel2abs(dw_ijlk[i], hcw_ijlk[i], dh_ijlk[i], WD_l[na, :, na])
+                if M == 1 or np.abs(dst_xyh_jlk[0][:, [0]] - dst_xyh_jlk[0]).max() > 1e-10:
+                    # destinations differ for differnet wd, use destinations calculated from
+                    # distances to include streamlines
+                    dst_xyh_jlk = _dst_xyh_jlk
+            if dst_xyh_jlk[0].shape[1:] == (1, 1) and M > J:
+                # same destination for all wd
+                x_j, y_j, h_j = [v[:, 0, 0] for v in dst_xyh_jlk]
+                lw_j, WS_eff_jlk, TI_eff_jlk = self.wfm._flow_map(x_j[:, na], y_j[:, na], h_j[:, na], sim_res.localWind,
+                                                                  wd_lst, ws_lst, sim_res)
+                _deficit_jlk = lw_j.WS_ilk - WS_eff_jlk
+
+                deficit_jmk = np.moveaxis([RGI([ws_lst], np.moveaxis(_deficit_jlk[:, l], 0, -1))(WS_eff_ilk[i, m])
+                                          for l, m in enumerate(m_lst)], -1, 0)
+            else:
+                assert dst_xyh_jlk[0].shape[2] == 1, "ExternalWFMWindFarm does not support inflow dependent positions"
+
+                def get_deficit_l(m):
+                    x_j, y_j, h_j = [np.broadcast_to(v_jlk, (J, L, K))[:, m] for v_jlk in dst_xyh_jlk]
+                    WD = WD_l[m]
+                    WS_eff = WS_eff_ilk[i, m]
+                    sr = sel_interp_all(sim_res)(dict(wd=WD, ws=WS_eff))
+                    lw_j, WS_eff_jlk, TI_eff_jlk = self.wfm._flow_map(x_j, y_j, h_j, sim_res.localWind, WD, WS_eff, sr)
+                    return lw_j.WS_ilk[:, 0] - WS_eff_jlk[:, 0]
+
+                deficit_jmk = np.moveaxis([get_deficit_l(m) for m in m_lst], 0, 1)
+            deficit_jlk = item_assign(deficit_jlk, m_lst, deficit_jmk, axis=1)
+        return deficit_jlk
+
+
+class ExternalXRAbsWindFarm(ExternalWindFarm):
+    def __init__(self, name, ds, windTurbines, wt_x, wt_y, relative_distance=True, include_wd=None):
+        ExternalWindFarm.__init__(self, name, windTurbines, wt_x, wt_y)
+        x = [ds.x.values, ds.y.values, ds.h.values, ds.wd.values, ds.ws.values]
+        self.deficit_interp = RGI(x, ds.deficit.values, bounds_error=False)
+        if include_wd is None:
+            include_wd = ds.wd.values
+        self.set_include_wd(include_wd)
+        self.relative_distance = relative_distance
+
+    def to_xarray(self):
+        di = self.deficit_interp
+        return xr.Dataset({'deficit': (['x', 'y', 'h', 'wd', 'ws'], di.values)},
+                          coords=dict(x=di.grid[0], y=di.grid[1], h=di.grid[2], wd=di.grid[3], ws=di.grid[4]))
+
+    @classmethod
+    def generate(cls, name, grid_xyh, windFarmModel, wt_x, wt_y, type=0, **kwargs):
+        sim_res = windFarmModel(wt_x, wt_y, type=type, **kwargs)
+        X, Y, H = np.meshgrid(*grid_xyh, indexing='ij')
+        x_j, y_j, h_j = X.flatten(), Y.flatten(), H.flatten()
+
+        lw_j, WS_eff_jlk, TI_eff_jlk = windFarmModel._flow_map(x_j[:, na], y_j[:, na], h_j[:, na], sim_res.localWind,
+                                                               sim_res.wd, sim_res.ws, sim_res)
+        deficit = lw_j.WS_ilk - WS_eff_jlk.reshape(X.shape + (WS_eff_jlk.shape[1:]))
+        ds = xr.Dataset({'deficit': (['x', 'y', 'h', 'wd', 'ws'], deficit)},
+                        coords=dict(x=grid_xyh[0], y=grid_xyh[1], h=np.atleast_1d(grid_xyh[2]), wd=sim_res.wd, ws=sim_res.ws))
+
+        from py_wake.site.streamline_distance import StreamlineDistance
+        return cls(name, ds, windFarmModel.windTurbines, wt_x, wt_y,
+                   relative_distance=not isinstance(windFarmModel.site.distance, StreamlineDistance))
+
+    def __call__(self, i, l, deficit_jlk, WS_eff_ilk, WS_ilk,
+                 dw_ijlk, hcw_ijlk, dh_ijlk,
+                 WD_ilk, IJLK, dst_xyh_jlk, **_):
+        I, J, L, K = IJLK
+        WD_l = WD_ilk[np.minimum(i, len(WD_ilk) - 1), :, 0]
+        m_lst = [m for m, wd in enumerate(WD_l) if l[m] and self.include_wd_func(wd)]
+        if deficit_jlk is None:
+            deficit_jlk = np.zeros((J, L, K))
+        if m_lst:
+            M = len(m_lst)
+            WS_eff_ilk[i]
+
+            if self.relative_distance:
+                # Simulation and flow map of external WF is computed with straight distance.
+                # The WF wakes must therefore be calcualted using the relative distances (including streamlines)
+                _dst_xyh_jlk = self.rel2abs(dw_ijlk[i], hcw_ijlk[i], dh_ijlk[i], WD_l[na, :, na])
+                if M == 1 or np.abs(dst_xyh_jlk[0][:, [0]] - dst_xyh_jlk[0]).max() > 1e-10:
+                    # destinations differ for different wd, use destinations calculated from
+                    # distances to include streamlines
+                    dst_xyh_jlk = _dst_xyh_jlk
+            if dst_xyh_jlk[0].shape[1:] == (1, 1) and M > J:
+                # same destination for all wd
+                # make interpolator with current destination and interpolate ws for each wd
+                di = self.deficit_interp
+                fm_jlk = RGI(di.grid[:3], di.values, bounds_error=False)(np.array([v[:, 0, 0] for v in dst_xyh_jlk]).T)
+                np.nan_to_num(fm_jlk, copy=False)
+                rgi_dst = RGI(di.grid[3:], np.moveaxis(fm_jlk, 0, -1), bounds_error=False)
+                WS_eff = WS_eff_ilk[i, m_lst]
+                WD = np.broadcast_to(WD_l[m_lst, na], WS_eff.shape)
+                deficit_jmk = np.moveaxis(
+                    rgi_dst(np.array([WD.flatten(), WS_eff.flatten()]).T), -1, 0).reshape((J, M, K))
+            else:
+                x = [np.broadcast_to(v[:, m_lst], (J, M, K)).flatten()
+                     for v in list(dst_xyh_jlk) + [WD_l[na, :, na], WS_eff_ilk[i, :][na]]]
+
+                deficit_jmk = self.deficit_interp(np.array(x).T).reshape((J, M, K))
+            np.nan_to_num(deficit_jmk, copy=False)
+            deficit_jlk = item_assign(deficit_jlk, m_lst, deficit_jmk, axis=1)
+        return deficit_jlk
+
+
+class ExternalXRRelWindFarm(ExternalXRAbsWindFarm):
+    @classmethod
+    def generate(cls, name, grid_xyh, windFarmModel, wt_x, wt_y, type=0, **kwargs):
+        sim_res = windFarmModel(wt_x, wt_y, type=type, **kwargs)
+        X, Y, H = np.meshgrid(*grid_xyh, indexing='ij')
+        dw, hcw, dh = X.flatten(), Y.flatten(), H.flatten()
+        wf_x = np.mean([np.min(wt_x), np.max(wt_x)])
+        wf_y = np.mean([np.min(wt_y), np.max(wt_y)])
+        wf_h = windFarmModel.windTurbines.hub_height(type)
+        deficit = []
+        for wd in tqdm(sim_res.wd.values, disable=1):
+            theta = np.deg2rad(270 - wd)
+            co, si = np.cos(theta), np.sin(theta)
+            x_j = co * dw - hcw * si + wf_x
+            y_j = si * dw + hcw * co + wf_y
+            h_j = dh + wf_h
+            lw_j, WS_eff_jlk, TI_eff_jlk = windFarmModel._flow_map(x_j[:, na], y_j[:, na], h_j[:, na], sim_res.localWind,
+                                                                   wd, sim_res.ws, sim_res)
+            deficit.append(lw_j.WS_ilk - WS_eff_jlk)
+        deficit = np.moveaxis(deficit, 0, 1).reshape(X.shape + (len(sim_res.wd), len(sim_res.ws)))
+        ds = xr.Dataset({'deficit': (['x', 'y', 'h', 'wd', 'ws'], deficit)},
+                        coords=dict(x=grid_xyh[0], y=grid_xyh[1], h=grid_xyh[2], wd=sim_res.wd, ws=sim_res.ws))
+
+        from py_wake.site.streamline_distance import StreamlineDistance
+        return cls(name, ds, windFarmModel.windTurbines, wt_x, wt_y,
+                   relative_distance=not isinstance(windFarmModel.site.distance, StreamlineDistance))
+
+    def __call__(self, i, l, deficit_jlk, WS_eff_ilk, WS_ilk,
+                 dw_ijlk, hcw_ijlk, dh_ijlk,
+                 WD_ilk, IJLK, dst_xyh_jlk, **_):
+        I, J, L, K = IJLK
+        WD_l = WD_ilk[np.minimum(i, len(WD_ilk) - 1), :, 0]
+        m_lst = [m for m, wd in enumerate(WD_l) if l[m] and self.include_wd_func(wd)]
+        if deficit_jlk is None:
+            deficit_jlk = np.zeros((J, L, K))
+        if not self.relative_distance:
+            xyh_ilk = [np.reshape(v, (1, 1, 1)) for v in [self.wf_x, self.wf_y, self.wf_h]]
+            dw_jlk, hcw_jlk, dh_jlk = [v[0] for v in StraightDistance()(*xyh_ilk, wd_l=WD_l, dst_xyh_jlk=dst_xyh_jlk)]
+        else:
+            dw_jlk, hcw_jlk, dh_jlk = dw_ijlk[i], hcw_ijlk[i], np.broadcast_to(dh_ijlk[i], (J, L, K))
+        if m_lst:
+            M = len(m_lst)
+            WS_eff_ilk[i]
+
+            x = [np.broadcast_to(v[:, m_lst], (J, M, K)).flatten()
+                 for v in [dw_jlk, hcw_jlk, dh_jlk,
+                           WD_l[na, :, na], WS_eff_ilk[i][na]]]
+
+            deficit_jmk = self.deficit_interp(np.array(x).T).reshape((J, M, K))
+            np.nan_to_num(deficit_jmk, copy=False)
+            deficit_jlk = item_assign(deficit_jlk, m_lst, deficit_jmk, axis=1)
+        return deficit_jlk

METADATA

 Metadata-Version: 2.4
 Name: py_wake
-Version: 2.6.13
+Version: 2.6.14
 Summary: Open source static wake modeling framework from DTU
@@ -5,0 +5,0 @@ Author: DTU Wind Energy

py_wake/deficit_models/deficit_model.py

@@ -115,3 +115,3 @@ from abc import ABC, abstractmethod
         R_ijlk = (D_src_il / 2)[:, na, :, na]
-        iw = ((dw_ijlk / R_ijlk >= -self.limiter) & (cabs(cw_ijlk) <= wake_radius_ijlk))
+        iw = ((dw_ijlk / R_ijlk >= -self.limiter) & (cabs(cw_ijlk) <= wake_radius_ijlk) & (wake_radius_ijlk > 0))
         deficit_ijlk = np.where(iw, 0., deficit_ijlk)
@@ -118,0 +118,0 @@

py_wake/flow_map.py

@@ -142,3 +142,3 @@ from py_wake import np
     def plot(self, data, clabel, levels=100, cmap=None, plot_colorbar=True, plot_windturbines=True,
-             normalize_with=1, ax=None):
+             normalize_with=1, ax=None, cax=None):
         """Plot data as contouf map
@@ -195,3 +195,3 @@
             if plot_colorbar:
-                plt.colorbar(c, label=clabel, ax=ax)
+                plt.colorbar(c, label=clabel, ax=ax, cax=cax)
         else:
@@ -203,2 +203,4 @@ raise NotImplementedError(
             self.plot_windturbines(normalize_with=normalize_with, ax=ax)
+            for ewf in self.windFarmModel.externalWindFarms:
+                ewf.plot(ax=ax)
         ax.axis('equal')
@@ -209,7 +211,8 @@ return c
         fm = self.windFarmModel
-
+        n_ewf = len(fm.externalWindFarms)
         # mean over wd and ws if present
-        x_i = self.simulationResult.x.mean(set(self.simulationResult.x.dims) - {'wt'}).values
-        y_i = self.simulationResult.y.mean(set(self.simulationResult.x.dims) - {'wt'}).values
-        type_i = self.simulationResult.type.data
+        n = len(self.simulationResult.wt) - n_ewf
+        x_i = self.simulationResult.x.mean(set(self.simulationResult.x.dims) - {'wt'}).values[:n]
+        y_i = self.simulationResult.y.mean(set(self.simulationResult.x.dims) - {'wt'}).values[:n]
+        type_i = self.simulationResult.type.values[:n]
         if self.plane[0] in ['XZ', "YZ"]:
@@ -256,3 +259,3 @@ h_i = self.simulationResult.h.values
     def plot_wake_map(self, levels=100, cmap=None, plot_colorbar=True, plot_windturbines=True,
-                      normalize_with=1, ax=None):
+                      normalize_with=1, ax=None, cax=None):
         """Plot effective wind speed contourf map
@@ -280,3 +283,3 @@
                          levels=levels, cmap=cmap, plot_colorbar=plot_colorbar,
-                         plot_windturbines=plot_windturbines, normalize_with=normalize_with, ax=ax)
+                         plot_windturbines=plot_windturbines, normalize_with=normalize_with, ax=ax, cax=cax)
 
@@ -512,2 +515,3 @@ def plot_ti_map(self, levels=100, cmap=None, plot_colorbar=True, plot_windturbines=True, ax=None):
     def __init__(self, x, y, h):
+        x, y, h = [np.atleast_1d(v) for v in [x, y, h]]
         assert len(x) == len(y) == len(h)
@@ -514,0 +518,0 @@ self.x = x

py_wake/noise_models/iso.py

@@ -56,3 +56,3 @@ import numpy as np
 
-    def ground_eff(self, ground_distance_ijlk, distance_ijlk, ground_type):
+    def ground_eff(self, ground_distance_ijlk, distance_ijlk, rec_h, ground_type):
         # Ground effects ISO
@@ -62,5 +62,7 @@
 
-        src_h_ijlk = self.distance.src_h_ilk[:, na]
-        rec_h_ijlk = self.distance.dst_h_j[na, :, na, na]
+        src_h_ilk = np.expand_dims(self.src_h, tuple(range(len(np.shape(self.src_h)), 3))) + .0
+        src_h_ijlk = src_h_ilk[:, na]
 
+        rec_h_ijlk = rec_h[na, :, na, na]
+
         hei_check_ijlk = 30.0 * (src_h_ijlk + rec_h_ijlk)
@@ -176,9 +178,9 @@
         rec_h = np.zeros_like(rec_x) + rec_h
-        rec_z = self.elevation_function(rec_x, rec_y)
-        self.distance.setup(self.src_x, self.src_y, self.src_h, self.src_z, [rec_x, rec_y, rec_h, rec_z])
-        ground_distance_ijlk = np.sqrt(self.distance.dx_ijlk**2 + self.distance.dy_ijlk**2)
-        distance_ijlk = np.sqrt(self.distance.dx_ijlk**2 + self.distance.dy_ijlk**2 + self.distance.dh_ijlk**2)
+        dx_ijlk, dy_ijlk, dh_ijlk = self.distance(self.src_x, self.src_y, self.src_h, wd_l=np.array([270]),
+                                                  dst_xyh_jlk=[rec_x, rec_y, rec_h])
+        ground_distance_ijlk = np.sqrt(dx_ijlk**2 + dy_ijlk**2)
+        distance_ijlk = np.sqrt(dx_ijlk**2 + dy_ijlk**2 + dh_ijlk**2)
 
         atm_abs_ijlkf = self.atmab(distance_ijlk, T0=Temp, RH0=RHum, p_s=patm)  # The atmospheric absorption term
-        ground_eff_ijlkf = self.ground_eff(ground_distance_ijlk, distance_ijlk, ground_type)
+        ground_eff_ijlkf = self.ground_eff(ground_distance_ijlk, distance_ijlk, rec_h, ground_type)
 
@@ -185,0 +187,0 @@ return ground_eff_ijlkf - atm_abs_ijlkf  # Delta_SPL

py_wake/rotor_avg_models/gaussian_overlap_model.py

@@ -10,2 +10,3 @@ import os
 from scipy.interpolate import RectBivariateSpline
+from py_wake.utils import gradients
 
@@ -25,3 +26,4 @@
             # dat = table.interp(R_sigma=R_sigma, CW_sigma=CW_sigma, method='cubic')
-            # performs a two-step interpolation which is slower and less accurate than the previous one-step interpolation
+            # performs a two-step interpolation which is slower and less accurate than
+            # the previous one-step interpolation
             dat = RectBivariateSpline(table.R_sigma.values, table.CW_sigma.values, table.values)(R_sigma, CW_sigma)
@@ -40,4 +42,5 @@ self._overlap_interpolator = GridInterpolator([R_sigma, CW_sigma], dat, bounds='limit')
         sigma_ijlk = func.__self__.sigma_ijlk(**kwargs)
-        overlap_factor_ijlk = self.overlap_interpolator(
-            np.array([((D_dst_ijl / 2)[..., na] / sigma_ijlk).flatten(), (cw_ijlk / sigma_ijlk).flatten()]).T).reshape(sigma_ijlk.shape)
+        r_sigma = ((D_dst_ijl / 2)[..., na] / sigma_ijlk).flatten()
+        cw_sigma = gradients.cabs(cw_ijlk / sigma_ijlk).flatten()
+        overlap_factor_ijlk = self.overlap_interpolator(np.array([r_sigma, cw_sigma]).T).reshape(sigma_ijlk.shape)
         return res_ijlk * overlap_factor_ijlk
@@ -44,0 +47,0 @@

py_wake/site/_site.py

@@ -214,5 +214,2 @@ import matplotlib.pyplot as plt
 
-    def wt2wt_distances(self, WD_ilk=None, wd_l=None):
-        return self.distance(WD_ilk=WD_ilk, wd_l=wd_l)
-
     @abstractmethod
@@ -219,0 +216,0 @@ def elevation(self, x_i, y_i):

py_wake/site/distance.py

@@ -21,2 +21,3 @@ from py_wake import np
         """
+        assert wind_direction in ['wd', 'WD_i'], "'StraightDistance.wind_direction must be 'wd' or 'WD_i'"
         self.wind_direction = wind_direction
@@ -30,6 +31,73 @@
 
-    def plot(self, WD_ilk=None, wd_l=None, src_idx=slice(None), dst_idx=slice(None)):
+    def __getstate__(self):
+        return {k: v for k, v in self.__dict__.items()
+                if k not in {'src_x_ilk', 'src_y_ilk', 'src_h_ilk', 'dst_x_j', 'dst_y_j', 'dst_h_j',
+                             'dx_iilk', 'dy_iilk', 'dh_iilk', 'dx_ijlk', 'dy_ijlk', 'dh_ij', 'src_eq_dst'}}
+
+    def get_pos(self, src_x_ilk, src_y_ilk, src_h_ilk, wd_l=None, WD_ilk=None, dst_xyh_jlk=None):
+        # ensure 3d and
+        # +.0 ensures float or complex
+        src_x_ilk, src_y_ilk, src_h_ilk = [np.expand_dims(v, tuple(range(len(np.shape(v)), 3))) + .0
+                                           for v in [src_x_ilk, src_y_ilk, src_h_ilk]]
+
+        if dst_xyh_jlk is None:
+            dst_x_jlk, dst_y_jlk, dst_h_jlk = src_x_ilk, src_y_ilk, src_h_ilk
+        else:
+            assert len(dst_xyh_jlk) == 3
+            dst_x_jlk, dst_y_jlk, dst_h_jlk = [np.expand_dims(v, tuple(range(len(np.shape(v)), 3))) + .0
+                                               for v in dst_xyh_jlk]
+        return (src_x_ilk, src_y_ilk, src_h_ilk), (dst_x_jlk, dst_y_jlk, dst_h_jlk)
+
+    def __call__(self, src_x_ilk, src_y_ilk, src_h_ilk, wd_l=None, WD_ilk=None, time=None, dst_xyh_jlk=None):
+        (src_x_ilk, src_y_ilk, src_h_ilk), (dst_x_jlk, dst_y_jlk, dst_h_jlk) = self.get_pos(
+            src_x_ilk, src_y_ilk, src_h_ilk, wd_l, WD_ilk, dst_xyh_jlk)
+        wd_l = np.asarray(wd_l)
+        if self.wind_direction == 'wd':
+            assert wd_l is not None, "wd_l must be specified when Distance.wind_direction='wd'"
+            WD_ilk = np.asarray(wd_l)[na, :, na]
+        else:
+            assert WD_ilk is not None, "WD_ilk must be specified when Distance.wind_direction='WD_i'"
+
+        cos_ilk, sin_ilk = self._cos_sin(WD_ilk)
+        if self.wind_direction == 'wd':
+
+            src_dw_ilk = -cos_ilk * src_x_ilk - sin_ilk * src_y_ilk
+            src_hcw_ilk = sin_ilk * src_x_ilk - cos_ilk * src_y_ilk
+
+            if dst_xyh_jlk is None:
+                dst_dw_jlk, dst_hcw_jlk, dst_h_jlk = src_dw_ilk, src_hcw_ilk, src_h_ilk
+            else:
+                cos_jlk, sin_jlk = self._cos_sin(wd_l[na, :, na])
+                dst_dw_jlk = (-cos_jlk * dst_x_jlk - sin_jlk * dst_y_jlk)
+                dst_hcw_jlk = (sin_jlk * dst_x_jlk - cos_jlk * dst_y_jlk)
+            dw_ijlk = dst_dw_jlk[na] - src_dw_ilk[:, na]
+            hcw_ijlk = dst_hcw_jlk[na] - src_hcw_ilk[:, na]
+        else:
+            dx_ijlk = dst_x_jlk[na] - src_x_ilk[:, na]
+            dy_ijlk = dst_y_jlk[na] - src_y_ilk[:, na]
+
+            cos_ilk, sin_ilk = self._cos_sin(WD_ilk)
+            dw_ijlk = -cos_ilk[:, na] * dx_ijlk - sin_ilk[:, na] * dy_ijlk
+            hcw_ijlk = sin_ilk[:, na] * dx_ijlk - cos_ilk[:, na] * dy_ijlk
+
+        dh_ijlk = dst_h_jlk[na, :] - src_h_ilk[:, na]
+
+        return dw_ijlk, hcw_ijlk, np.broadcast_to(dh_ijlk, hcw_ijlk.shape)
+
+    def dw_order_indices(self, src_x_ilk, src_y_ilk, wd_l):
+        WD_ilk = np.asarray(wd_l)[na, :, na]
+        src_x_ilk, src_y_ilk = [np.expand_dims(v, tuple(range(len(np.shape(v)), 3))) + .0
+                                for v in [src_x_ilk, src_y_ilk]]
+
+        cos_ilk, sin_ilk = self._cos_sin(WD_ilk)
+
+        src_dw_ilk = -cos_ilk * src_x_ilk - sin_ilk * src_y_ilk
+        return np.moveaxis(np.argsort(src_dw_ilk, 0), 0, -1)
+
+    def plot(self, src_x_ilk, src_y_ilk, src_h_ilk, wd_l=None, WD_ilk=None, dst_xyh_jlk=None):
         import matplotlib.pyplot as plt
+        (src_x_ilk, src_y_ilk, src_h_ilk), (dst_x_jlk, dst_y_jlk, dst_h_jlk) = self.get_pos(
+            src_x_ilk, src_y_ilk, src_h_ilk, wd_l, WD_ilk, dst_xyh_jlk)
 
-        dw_ijlk, hcw_ijlk, _ = self(WD_ilk=WD_ilk, wd_l=wd_l)
+        dw_ijlk, hcw_ijlk, _ = self(src_x_ilk, src_y_ilk, src_h_ilk, WD_ilk=WD_ilk, wd_l=wd_l, dst_xyh_jlk=dst_xyh_jlk)
         if self.wind_direction == 'wd':
@@ -47,8 +115,7 @@ WD_ilk = np.asarray(wd_l)[na, :, na]
             f = 2
-            for i, x_, y_ in zip(np.arange(len(self.src_x_ilk))[
-                                 src_idx], self.src_x_ilk[src_idx, 0, 0], self.src_y_ilk[src_idx, 0, 0]):
+            for i, x_, y_ in zip(np.arange(len(src_x_ilk)), src_x_ilk[:, 0, 0], src_y_ilk[:, 0, 0]):
                 c = colors[i % len(colors)]
                 ax.plot(x_, y_, '2', color=c, ms=10, mew=3)
-                for j, dst_x, dst_y in zip(np.arange(len(self.dst_x_j))[dst_idx],
-                                           self.dst_x_j[dst_idx, 0, 0], self.dst_y_j[dst_idx, 0, 0]):
+                for j, dst_x, dst_y in zip(np.arange(len(dst_x_jlk)),
+                                           dst_x_jlk[:, 0, 0], dst_y_jlk[:, 0, 0]):
                     ax.arrow(x_ - j / f, y_ - j / f, -np.cos(theta) * dw_ijlk[i, j, l, 0], -
@@ -58,77 +125,7 @@ np.sin(theta) * dw_ijlk[i, j, l, 0], width=.3, color=c)
                             [dst_y - i / f, dst_y + np.cos(theta) * hcw_ijlk[i, j, l, 0] - i / f], '--', color=c)
-            plt.plot(self.src_x_ilk[:, 0, 0], self.src_y_ilk[:, 0, 0], 'k2')
+            plt.plot(src_x_ilk[:, 0, 0], src_y_ilk[:, 0, 0], 'k2')
+            plt.plot(dst_x_jlk[:, 0, 0], dst_y_jlk[:, 0, 0], 'k2')
             ax.axis('equal')
 
-    def setup(self, src_x_ilk, src_y_ilk, src_h_ilk, src_z_ilk, dst_xyhz_j=None):
-        # ensure 3d and
-        # +.0 ensures float or complex
-        src_x_ilk, src_y_ilk, src_h_ilk, src_z_ilk = [np.expand_dims(v, tuple(range(len(np.shape(v)), 3))) + .0
-                                                      for v in [src_x_ilk, src_y_ilk, src_h_ilk, src_z_ilk]]
-        self.src_x_ilk, self.src_y_ilk, self.src_h_ilk = src_x_ilk, src_y_ilk, src_h_ilk
 
-        self.dx_iilk = src_x_ilk - src_x_ilk[:, na]
-        self.dy_iilk = src_y_ilk - src_y_ilk[:, na]
-        self.dh_iilk = src_h_ilk - src_h_ilk[:, na]
-        self.dz_iilk = src_z_ilk - src_z_ilk[:, na]
-        if dst_xyhz_j is None:
-            dst_x_j, dst_y_j, dst_h_j, dst_z_j = src_x_ilk, src_y_ilk, src_h_ilk, src_z_ilk
-            self.dx_ijlk, self.dy_ijlk, self.dh_ijlk, self.dz_ijlk = self.dx_iilk, self.dy_iilk, self.dh_iilk, self.dz_iilk
-            self.src_eq_dst = True
-        else:
-            dst_x_j, dst_y_j, dst_h_j, dst_z_j = map(np.asarray, dst_xyhz_j)
-            self.dx_ijlk = dst_x_j[na, :, na, na] - src_x_ilk[:, na]
-            self.dy_ijlk = dst_y_j[na, :, na, na] - src_y_ilk[:, na]
-            self.dh_ijlk = dst_h_j[na, :, na, na] - src_h_ilk[:, na]
-            self.dz_ijlk = dst_z_j[na, :, na, na] - src_z_ilk[:, na]
-            self.src_eq_dst = False
-        self.dst_x_j, self.dst_y_j, self.dst_h_j, self.dst_z_j = dst_x_j, dst_y_j, dst_h_j, dst_z_j
-
-    def __getstate__(self):
-        return {k: v for k, v in self.__dict__.items()
-                if k not in {'src_x_ilk', 'src_y_ilk', 'src_h_ilk', 'dst_x_j', 'dst_y_j', 'dst_h_j',
-                             'dx_iilk', 'dy_iilk', 'dh_iilk', 'dx_ijlk', 'dy_ijlk', 'dh_ij', 'src_eq_dst'}}
-
-    def __call__(self, WD_ilk=None, wd_l=None, src_idx=slice(None), dst_idx=slice(None)):
-        assert hasattr(self, 'dx_ijlk'), "wind_direction setup must be called first"
-        assert self.wind_direction in ['wd', 'WD_i'], "'StraightDistance.wind_direction must be 'wd' or 'WD_i'"
-
-        if self.wind_direction == 'wd':
-            assert wd_l is not None, "wd_l must be specified when Distance.wind_direction='wd'"
-            WD_ilk = np.asarray(wd_l)[na, :, na]
-        else:
-            assert WD_ilk is not None, "WD_ilk must be specified when Distance.wind_direction='WD_i'"
-
-        if len(np.shape(dst_idx)) == 2:
-            # dst_idx depends on wind direction
-            cos_jlk, sin_jlk = self._cos_sin(WD_ilk[0, na])
-            i_wd_l = np.arange(np.shape(dst_idx)[1])
-            dx_jlk = self.dx_iilk[src_idx, dst_idx, np.minimum(i_wd_l, self.dx_iilk.shape[2] - 1).astype(int)]
-            dy_jlk = self.dy_iilk[src_idx, dst_idx, np.minimum(i_wd_l, self.dy_iilk.shape[2] - 1).astype(int)]
-            dh_jlk = self.dh_iilk[src_idx, dst_idx, np.minimum(i_wd_l, self.dh_iilk.shape[2] - 1).astype(int)]
-            dw_jlk = (-cos_jlk * dx_jlk - sin_jlk * dy_jlk)
-            hcw_jlk = (sin_jlk * dx_jlk - cos_jlk * dy_jlk)
-            return dw_jlk[na], hcw_jlk[na], dh_jlk[na]
-        else:
-            # dst_idx independent of wind direction
-            cos_ijlk, sin_ijlk = self._cos_sin(WD_ilk[:, na])
-            dx_ijlk = self.dx_ijlk[src_idx][:, dst_idx]
-            dy_ijlk = self.dy_ijlk[src_idx][:, dst_idx]
-
-            dw_ijlk = -cos_ijlk * dx_ijlk - sin_ijlk * dy_ijlk
-            hcw_ijlk = sin_ijlk * dx_ijlk - cos_ijlk * dy_ijlk
-            # +0 ~ autograd safe copy (broadcast_to returns readonly array)
-            dh_ijlk = np.broadcast_to(self.dh_ijlk[src_idx][:, dst_idx], dw_ijlk.shape) + 0
-            return dw_ijlk, hcw_ijlk, dh_ijlk
-
-    def dw_order_indices(self, wd_l):
-        assert hasattr(self, 'dx_ijlk'), "method setup must be called first"
-        I, J, *_ = self.dx_ijlk.shape
-        assert I == J
-        cos_l, sin_l = self._cos_sin(np.asarray(wd_l))
-        # return np.argsort(-cos_l[:, na] * np.asarray(src_x_ilk)[na] - sin_l[:, na] * np.asarray(src_y_ilk)[na], 1)
-        dw_iil = -cos_l[na, na, :, na] * self.dx_ijlk - sin_l[na, na, :, na] * self.dy_ijlk
-        dw_order_indices_lkd = np.moveaxis(np.argsort((dw_iil > 0).sum(0), 0), 0, -1)
-        return dw_order_indices_lkd
-
-
 class TerrainFollowingDistance(StraightDistance):
@@ -139,34 +136,27 @@ def __init__(self, distance_resolution=1000, wind_direction='wd', **kwargs):
 
-    def setup(self, src_x_ilk, src_y_ilk, src_h_ilk, src_z_ilk, dst_xyhz_j=None):
-        StraightDistance.setup(self, src_x_ilk, src_y_ilk, src_h_ilk, src_z_ilk, dst_xyhz_j=dst_xyhz_j)
-        if len(src_x_ilk) == 0:
-            return
+    def __call__(self, src_x_ilk, src_y_ilk, src_h_ilk,
+                 WD_ilk=None, wd_l=None, time=None, dst_xyh_jlk=None):
+        # StraightDistance.setup(self, src_x_ilk, src_y_ilk, src_h_ilk, src_z_ilk, dst_xyh_jlk=dst_xyh_jlk)
+        # if len(src_x_ilk) == 0:
+        #     return
         # Calculate distance between src and dst and project to the down wind direction
-        assert self.src_x_ilk.shape[1:] == (
-            1, 1), 'TerrainFollowingDistance does not support flowcase dependent positions'
-        src_x_ilk, src_y_ilk, src_h_ilk = self.src_x_ilk[:, 0, 0], self.src_y_ilk[:, 0, 0], self.src_h_ilk[:, 0, 0]
-        if len(self.dst_x_j.shape) == 3:
-            dst_x_j, dst_y_j = self.dst_x_j[:, 0, 0], self.dst_y_j[:, 0, 0]
-        else:
-            dst_x_j, dst_y_j = self.dst_x_j, self.dst_y_j
 
+        (src_x_ilk, src_y_ilk, src_h_ilk), (dst_x_jlk, dst_y_jlk, dst_h_jlk) = self.get_pos(
+            src_x_ilk, src_y_ilk, src_h_ilk, wd_l, WD_ilk, dst_xyh_jlk)
+        dw_ijlk, hcw_ijlk, dh_ijlk = StraightDistance.__call__(self, src_x_ilk, src_y_ilk, src_h_ilk,
+                                                               WD_ilk=WD_ilk, wd_l=wd_l, dst_xyh_jlk=dst_xyh_jlk)
+
+        assert src_x_ilk.shape[2] == 1, 'TerrainFollowingDistance does not support flowcase dependent positions'
+
+        src_x_i, src_y_i = src_x_ilk[:, 0, 0], src_y_ilk[:, 0, 0]
+
+        dst_x_j, dst_y_j = dst_x_jlk[:, 0, 0], dst_y_jlk[:, 0, 0]
+
         # Generate interpolation lines
-
-        if (self.src_eq_dst and self.dx_ijlk.shape[0] > 1):
-            # calculate upper triangle of d_ij(distance from i to j) only
-            xy = np.array([(np.linspace(src_x, dst_x, self.distance_resolution),
-                            np.linspace(src_y, dst_y, self.distance_resolution))
-                           for i, (src_x, src_y) in enumerate(zip(src_x_ilk, src_y_ilk))
-                           for dst_x, dst_y in zip(dst_x_j[i + 1:], dst_y_j[i + 1:])])
-            upper_tri_only = True
-            self.theta_ij = gradients.arctan2(dst_y_j[na, :] - src_y_ilk[:, na],
-                                              dst_x_j[na, :] - src_x_ilk[:, na])
-        else:
-            xy = np.array([(np.linspace(src_x, dst_x, self.distance_resolution),
-                            np.linspace(src_y, dst_y, self.distance_resolution))
-                           for src_x, src_y in zip(src_x_ilk, src_y_ilk)
-                           for dst_x, dst_y in zip(dst_x_j, dst_y_j)])
-            self.theta_ij = gradients.arctan2(dst_y_j[na, :, ] - src_y_ilk[:, na],
-                                              dst_x_j[na, :] - src_x_ilk[:, na])
-            upper_tri_only = False
+        xy = np.array([(np.linspace(src_x, dst_x, self.distance_resolution),
+                        np.linspace(src_y, dst_y, self.distance_resolution))
+                       for src_x, src_y in zip(src_x_i, src_y_i)
+                       for dst_x, dst_y in zip(dst_x_j, dst_y_j)])
+        theta_ij = gradients.arctan2(dst_y_j[na, :, ] - src_y_i[:, na],
+                                     dst_x_j[na, :] - src_x_i[:, na])
         x, y = xy[:, 0], xy[:, 1]
@@ -182,17 +172,4 @@
 
-        if upper_tri_only:
-            # d_ij = np.zeros(self.dx_ijlk.shape)
-            # d_ij[np.triu(np.eye(len(src_x_ilk)) == 0)] = s  # set upper triangle
+        d_ij = s.reshape(dw_ijlk.shape[:2])
 
-            # same as above without item assignment
-            n = len(src_x_ilk)
-            d_ij = np.array([np.concatenate([[0] * (i + 1),
-                                             s[int(n * i - (i * (i + 1) / 2)):][:n - i - 1]])
-                             for i in range(n)])  # set upper and lower triangle
-            d_ij += d_ij.T
-        else:
-            d_ij = s.reshape(self.dx_ijlk.shape[:2])
-        self.d_ij = d_ij
-
-    def __call__(self, WD_ilk=None, wd_l=None, src_idx=slice(None), dst_idx=slice(None)):
         # project terrain following distance between wts onto downwind direction
@@ -202,4 +179,2 @@ # instead of projecting the distances onto first x,y and then onto down wind direction
 
-        _, hcw_ijlk, dh_ijlk = StraightDistance.__call__(self, WD_ilk=WD_ilk, wd_l=wd_l,
-                                                         src_idx=src_idx, dst_idx=dst_idx)
         if self.wind_direction == 'wd':
@@ -209,17 +184,7 @@ WD_ilk = np.asarray(wd_l)[na, :, na]
         WD_il = mean_deg(WD_ilk, 2)
+        dir_ij = 90 - rad2deg(theta_ij)
+        wdir_offset_ijl = np.asarray(WD_il)[:, na] - dir_ij[:, :, na]
+        theta_ijl = deg2rad(90 - wdir_offset_ijl)
+        dw_ijlk = (- np.sin(theta_ijl) * d_ij[:, :, na])[..., na]
 
-        if len(np.shape(dst_idx)) == 2:
-            # dst_idx depends on wind direction
-            WD_l = WD_il[0]
-            dir_jl = 90 - rad2deg(self.theta_ij[src_idx, dst_idx])
-            wdir_offset_jl = WD_l[na, :] - dir_jl
-            theta_jl = deg2rad(90 - wdir_offset_jl)
-
-            dw_ijlk = (- np.sin(theta_jl) * self.d_ij[src_idx, dst_idx])[na, :, :, na]
-        else:
-            dir_ij = 90 - rad2deg(self.theta_ij[src_idx, ][:, dst_idx])
-            wdir_offset_ijl = np.asarray(WD_il)[:, na] - dir_ij[:, :, na]
-            theta_ijl = deg2rad(90 - wdir_offset_ijl)
-            dw_ijlk = (- np.sin(theta_ijl) * self.d_ij[src_idx][:, dst_idx][:, :, na])[..., na]
-
         return dw_ijlk, hcw_ijlk, dh_ijlk

py_wake/site/jit_streamline_distance.py

@@ -1,108 +0,3 @@
-from py_wake.site.distance import StraightDistance
-from py_wake import np
-from numpy import newaxis as na
-import matplotlib.pyplot as plt
-from py_wake.examples.data.ParqueFicticio._parque_ficticio import ParqueFicticioSite
-from py_wake.examples.data.hornsrev1 import V80
-from py_wake.deficit_models.noj import NOJ
-from py_wake.flow_map import XYGrid
-from py_wake.utils.streamline import VectorField3D
+import warnings
 
-
-class JITStreamlineDistance(StraightDistance):
-    """Just-In-Time Streamline Distance
-    Calculates downwind crosswind and vertical distance along streamlines.
-    Streamlines calculated in each call
-    """
-
-    def __init__(self, vectorField, step_size=20):
-        """Parameters
-        ----------
-        vectorField : VectorField3d
-        step_size : int for float
-            Size of linear streamline steps
-        """
-        StraightDistance.__init__(self, wind_direction='wd')
-        self.vectorField = vectorField
-        self.step_size = step_size
-
-    def __call__(self, wd_l, WD_ilk, src_idx=slice(None), dst_idx=slice(None)):
-        assert self.src_x_ilk.shape[1:] == (1, 1), 'JITStreamlineDistance does not support flowcase dependent positions'
-        start_points_m = np.array([self.src_x_ilk[src_idx], self.src_y_ilk[src_idx],
-                                   self.src_h_ilk[src_idx]])[:, :, 0, 0].T
-
-        if len(np.shape(dst_idx)) == 2:
-            # dst_idx depends on wind direction
-            dw_jl, hcw_jl, dh_jl = [v[0, :, :, 0] for v in StraightDistance.__call__(self, wd_l=wd_l,
-                                                                                     src_idx=src_idx, dst_idx=dst_idx)]
-            dw_mj, hcw_mj, dh_mj = [np.moveaxis(v, 0, 1) for v in [dw_jl, hcw_jl, dh_jl]]
-
-            i_wd_l = np.arange(np.shape(dst_idx)[1])
-            dz_jlk = self.dz_iilk[src_idx, dst_idx, np.minimum(i_wd_l, self.dh_iilk.shape[2] - 1).astype(int)]
-            dz_mj = np.moveaxis(dz_jlk[:, :, 0], 0, 1)
-
-            wd_m = wd_l
-        else:
-            # dst_idx independent of wind direction
-            dw_ijlk, hcw_ijlk, dh_ijlk = StraightDistance.__call__(self, wd_l=wd_l,
-                                                                   src_idx=src_idx, dst_idx=dst_idx)
-            # +0 ~ autograd safe copy (broadcast_to returns readonly array)
-            dz_ijlk = np.broadcast_to(self.dz_ijlk[src_idx][:, dst_idx], dw_ijlk.shape) + 0
-            I, J, L, K = dw_ijlk.shape
-            dw_mj, hcw_mj, dh_mj, dz_mj = [np.moveaxis(v, 1, 2).reshape(I * L, J)
-                                           for v in [dw_ijlk, hcw_ijlk, dh_ijlk, dz_ijlk]]
-
-            wd_m = np.tile(wd_l, I)
-            start_points_m = np.repeat(start_points_m, L, 0)
-
-        stream_lines = self.vectorField.stream_lines(wd_m, start_points=start_points_m, dw_stop=dw_mj.max(1),
-                                                     step_size=self.step_size)
-
-        dxyz = np.diff(np.concatenate([stream_lines[:, :1], stream_lines], 1), 1, -2)
-        length_is = np.cumsum(np.sqrt(np.sum(dxyz**2, -1)), -1)
-        dist_xyz = stream_lines - start_points_m[:, na]
-        t = np.deg2rad(270 - wd_m)[:, na]
-        dw_is = dist_xyz[:, :, 0] * np.cos(t) + dist_xyz[:, :, 1] * np.sin(t)
-        hcw_is = dist_xyz[:, :, 0] * np.sin(t) - dist_xyz[:, :, 1] * np.cos(t)
-
-        for m, (dw_j, dw_s, hcw_s, dh_s, length_s) in enumerate(
-                zip(dw_mj, dw_is, hcw_is, dist_xyz[:, :, 2], length_is)):
-            dw = dw_j > 0
-            hcw_mj[m, dw] += np.interp(dw_j[dw], dw_s, hcw_s)
-            dh_mj[m, dw] -= np.interp(dw_j[dw], dw_s, dh_s)
-            dw_mj[m, dw] = np.interp(dw_j[dw], dw_s, length_s)
-
-        # streamline dh contains absolute height different, but pywake needs differences relative to ground, so
-        # we need to subtract elevation differences, dz
-        dh_mj += dz_mj
-
-        if len(np.shape(dst_idx)) == 2:
-
-            return np.array([np.moveaxis(v, 0, 1)[na, :, :, na] for v in [dw_mj, hcw_mj, dh_mj]])
-        else:
-            return [v.reshape((I, J, L, K)) for v in [dw_mj, hcw_mj, dh_mj]]
-
-
-def main():
-    if __name__ == '__main__':
-
-        wt = V80()
-        vf3d = VectorField3D.from_WaspGridSite(ParqueFicticioSite())
-        site = ParqueFicticioSite(distance=JITStreamlineDistance(vf3d))
-
-        x, y = site.initial_position[:].T
-        wfm = NOJ(site, wt)
-        wd = 330
-        sim_res = wfm(x, y, wd=[wd], ws=10)
-        fm = sim_res.flow_map(XYGrid(x=np.linspace(site.ds.x[0].item(), site.ds.x[-1].item(), 500),
-                                     y=np.linspace(site.ds.y[0].item(), site.ds.y[-1].item(), 500)))
-        stream_lines = vf3d.stream_lines(wd=np.full(x.shape, wd), start_points=np.array([x, y, np.full(x.shape, 70)]).T,
-                                         dw_stop=y - 6504700)
-        fm.plot_wake_map()
-        for sl in stream_lines:
-            plt.plot(sl[:, 0], sl[:, 1])
-
-        plt.show()
-
-
-main()
+warnings.warn('jit_streamline_distance.JITStreamLineDistance has been renamed to streamline_distance.StreamLineDistance', DeprecationWarning)

py_wake/tests/test_deficit_models/test_deficit_models.py

@@ -31,3 +31,3 @@ import pytest
 from py_wake.deficit_models.utils import ct2a_mom1d
-from py_wake.wind_turbines._wind_turbines import WindTurbine, WindTurbines
+from py_wake.wind_turbines._wind_turbines import WindTurbines
 from py_wake.site._site import UniformSite
@@ -67,6 +67,6 @@ from py_wake.deficit_models.rans_lut import RANSLUTDemoDeficit
                      smooth2zero_x=250, smooth2zero_y=50),
-         (404307.913694007, [9962.5736, 9801.78637, 12608.917, 15452.89633, 22577.55462,
-                             27901.06282, 43479.02414, 49825.74736, 25105.68548, 15542.67624,
-                             16918.79822, 35640.98079, 76856.89565, 19752.83273, 13882.09085,
-                             8998.39151])),
+         (404307.96007812454, [9962.59462, 9801.78637, 12608.917, 15452.89633, 22577.55462,
+                               27901.06281, 43479.02414, 49825.74736, 25105.71085, 15542.67624,
+                               16918.79822, 35640.98079, 76856.89565, 19752.83273, 13882.09085,
+                               8998.39151])),
         (GCLDeficit(), (370863.6246093183,
@@ -232,4 +232,4 @@ [9385.75387, 8768.52105, 11450.13309, 14262.42186, 21178.74926,
     y_j = np.linspace(-1500, 1500, 100)
-
-    flow_map = wf_model(x, y, wd=0, ws=9).flow_map(HorizontalGrid(x_j, y_j))
+    wf_model(x, y, wd=1e-6, ws=9).flow_map(HorizontalGrid(x_j[100:133:2], y_j[[49]]))
+    flow_map = wf_model(x, y, wd=1e-6, ws=9).flow_map(HorizontalGrid(x_j, y_j))
     X, Y = flow_map.X, flow_map.Y
@@ -349,17 +349,17 @@ Z = flow_map.WS_eff_xylk[:, :, 0, 0]
     # test that the result is equal to last run (no evidens that  these number are correct)
-    [(BastankhahGaussianDeficit(ct2a=ct2a_mom1d), (345846.3355259293,
-                                                   [8835.30563, 7877.90062, 11079.66832, 13565.65235, 18902.99769,
-                                                    24493.53897, 38205.75284, 40045.9948, 22264.97018, 12662.90784,
-                                                    14650.96535, 31289.90349, 65276.92307, 17341.39229, 12021.3049,
-                                                    7331.15717])),
+    [(BastankhahGaussianDeficit(ct2a=ct2a_mom1d), (345844.56385386083,
+                                                   [8835.30376, 7877.89335, 11079.53783, 13565.45008, 18902.94686,
+                                                    24493.17471, 38205.304, 40045.95482, 22264.96545, 12662.88276,
+                                                    14650.96516, 31289.89954, 65276.44783, 17341.3905, 12021.30506,
+                                                    7331.14214])),
      (ZongGaussianDeficit(ct2a=ct2a_mom1d, eps_coeff=0.35),
-      (342944.4057168523, [8674.44232, 7806.22033, 11114.78804, 13549.48197, 18895.50866,
-                           24464.34244, 38326.85532, 39681.61999, 21859.59465, 12590.25899,
-                           14530.24656, 31158.81189, 63812.14454, 17268.73912, 11922.25359,
-                           7289.09731])),
+      (342945.18621135753, [8674.44233, 7806.22057, 11114.78784, 13549.48164, 18895.50866,
+                            24464.34303, 38326.85602, 39681.61877, 21860.37542, 12590.25945,
+                            14530.24631, 31158.81143, 63812.14454, 17268.73937, 11922.2538,
+                            7289.09705])),
      (NiayifarGaussianDeficit(ct2a=ct2a_mom1d),
-      (349044.6891842835, [8888.36909, 8025.38191, 11134.3202, 13643.08009, 19503.25093,
-                           24633.33906, 38394.20758, 40795.69138, 22398.69011, 12972.04355,
-                           14504.45911, 31328.69308, 66049.04782, 17362.89014, 11901.09465,
-                           7510.13048])),
+      (349044.6889144409, [8888.3691, 8025.38213, 11134.32002, 13643.07981, 19503.25093,
+                           24633.33957, 38394.20821, 40795.69026, 22398.6901, 12972.04395,
+                           14504.45888, 31328.69272, 66049.04782, 17362.89034, 11901.09484,
+                           7510.13025])),
 
@@ -375,3 +375,3 @@ ])
 
-    aep_ilk = wf_model(x, y, wd=np.arange(0, 360, 22.5), ws=[9.8]).aep_ilk(normalize_probabilities=True)
+    aep_ilk = wf_model(x, y, wd=np.arange(0, 360, 22.5) + 1e-6, ws=[9.8]).aep_ilk(normalize_probabilities=True)
     aep_MW_l = aep_ilk.sum((0, 2)) * 1000
@@ -495,4 +495,4 @@
                                    *args, **kwargs),
-      (404411.56948244764, [9963.95691, 9804.76072, 12613.39798, 15457.51387, 22582.08811,
-                            27909.40005, 43494.4758, 49840.86701, 25109.17142, 15546.80155,
+      (404411.61587925453, [9963.97791, 9804.76072, 12613.39798, 15457.51387, 22582.08811,
+                            27909.40005, 43494.4758, 49840.86701, 25109.19683, 15546.80155,
                             16923.11586, 35647.55255, 76875.57937, 19756.4749, 13885.63352,
@@ -499,0 +499,0 @@ 9000.77984])),

py_wake/tests/test_deficit_models/test_rans_lut.py

@@ -37,10 +37,11 @@ from py_wake import np
     y_j = np.linspace(-1500, 1500, 300)
-    flow_map70 = simres.flow_map(HorizontalGrid(x_j, y_j, h=70))
-    flow_map73 = simres.flow_map(HorizontalGrid(x_j, y_j, h=73))
 
-    X, Y = flow_map70.XY
-    Z70 = flow_map70.WS_eff_xylk[:, :, 0, 0]
-    Z73 = flow_map73.WS_eff_xylk[:, :, 0, 0]
+    if 0:
+        flow_map70 = simres.flow_map(HorizontalGrid(x_j, y_j, h=70))
+        flow_map73 = simres.flow_map(HorizontalGrid(x_j, y_j, h=73))
 
-    if 0:
+        X, Y = flow_map70.XY
+        Z70 = flow_map70.WS_eff_xylk[:, :, 0, 0]
+        Z73 = flow_map73.WS_eff_xylk[:, :, 0, 0]
+
         flow_map70.plot_wake_map(levels=np.arange(6, 10.5, .1))
@@ -52,8 +53,27 @@ plt.plot(X[0], Y[140])
         plt.plot(X[0, 100:400:10], Z70[140, 100:400:10], '.')
-        print(list(np.round(Z70.data[140, 100:400:10], 4)))
-        print(list(np.round(Z73.data[140, 100:400:10], 4)))
+        print(np.round(Z70.data[140, 100:400:10], 4).tolist())
+        print(np.round(Z73.data[140, 100:400:10], 4).tolist())
         plt.legend()
         plt.show()
 
+    flow_map70 = simres.flow_map(HorizontalGrid(x_j[100:400:10], y_j[[140]], h=70))
+    flow_map73 = simres.flow_map(HorizontalGrid(x_j[100:400:10], y_j[[140]], h=73))
 
+    X, Y = flow_map70.XY
+    Z70 = flow_map70.WS_eff_xylk[0, :, 0, 0]
+    Z73 = flow_map73.WS_eff_xylk[0, :, 0, 0]
+
+    npt.assert_array_almost_equal(
+        Z70,
+        [10.0199, 10.0235, 10.0293, 10.0358, 9.8078, 5.5094, 4.718, 9.4776, 10.034, 10.0113, 10.0116, 10.027, 10.0617,
+         9.4429, 5.6, 8.8807, 10.0984, 10.0073, 9.9856, 9.96, 9.0589, 7.4311, 3.7803, 6.5824, 10.1399, 10.0271,
+         9.9978, 9.9919, 9.9907, 9.9906], 4)
+
+    npt.assert_array_almost_equal(
+        Z73,
+        [10.0198, 10.0234, 10.0291, 10.0351, 9.8118, 5.6347, 4.7181, 9.4902, 10.0332, 10.0113, 10.0116, 10.0267,
+         10.0607, 9.4732, 5.5569, 8.9562, 10.0969, 10.0073, 9.9858, 9.9601, 9.0821, 7.4662, 3.8426, 6.6404, 10.1375,
+         10.0268, 9.9979, 9.992, 9.9908, 9.9906], 4)
+
+
 def test_rans_lut():
@@ -76,23 +96,3 @@ # move turbine 1 600 300
                                                             0.80567201, 0.80495685])
-    x_j = np.linspace(-1500, 1500, 500)
-    y_j = np.linspace(-1500, 1500, 300)
-    flow_map70 = simres.flow_map(HorizontalGrid(x_j, y_j, h=70))
-    flow_map73 = simres.flow_map(HorizontalGrid(x_j, y_j, h=73))
 
-    X, Y = flow_map70.XY
-    Z70 = flow_map70.WS_eff_xylk[:, :, 0, 0]
-    Z73 = flow_map73.WS_eff_xylk[:, :, 0, 0]
-
-    if 0:
-        flow_map70.plot_wake_map(levels=np.arange(6, 10.5, .1))
-        plt.plot(X[0], Y[140])
-        plt.figure()
-        plt.plot(X[0], Z70[140, :], label="Z=70m")
-        plt.plot(X[0], Z73[140, :], label="Z=73m")
-        plt.plot(X[0, 100:400:10], Z70[140, 100:400:10], '.')
-        print(list(np.round(Z70.data[140, 100:400:10], 4)))
-        print(list(np.round(Z73.data[140, 100:400:10], 4)))
-        plt.legend()
-        plt.show()
-
     # Test Power as RANS post step
@@ -157,3 +157,4 @@ dataset = xr.open_dataset(demo_lut)
     # Test MOST shear, stable.
-    # WS_eff_star does not change since we not change the actual simulation using Site with MOSTShear and LUTs with stability
+    # WS_eff_star does not change since we not change the actual simulation
+    # using Site with MOSTShear and LUTs with stability
     aDControl2 = ADControl.from_lut([dataset, lut2], wts, ws_cutin=4, ws_cutout=25, dws=1.0, cal_TI=0.06, cal_zeta=0.5)
@@ -215,23 +216,2 @@ ellipsys_power, WS_eff_star, ct_star = get_Ellipsys_equivalent_output(simres, aDControl2)
 
-    x_j = np.linspace(-1500, 1500, 500)
-    y_j = np.linspace(-1500, 1500, 300)
-    flow_map70 = simres.flow_map(HorizontalGrid(x_j, y_j, h=70), wd=30, ws=10)
-    flow_map73 = simres.flow_map(HorizontalGrid(x_j, y_j, h=73), wd=30, ws=10)
-
-    X, Y = flow_map70.XY
-    Z70 = flow_map70.WS_eff_xylk[:, :, 0, 0]
-    Z73 = flow_map73.WS_eff_xylk[:, :, 0, 0]
-
-    if 0:
-        flow_map70.plot_wake_map(levels=np.arange(6, 10.5, .1))
-        plt.plot(X[0], Y[140])
-        plt.figure()
-        plt.plot(X[0], Z70[140, :], label="Z=70m")
-        plt.plot(X[0], Z73[140, :], label="Z=73m")
-        plt.plot(X[0, 100:400:10], Z70[140, 100:400:10], '.')
-        print(list(np.round(Z70.data[140, 100:400:10], 4)))
-        print(list(np.round(Z73.data[140, 100:400:10], 4)))
-        plt.legend()
-        plt.show()
-
     # Test Power as RANS post step
@@ -238,0 +218,0 @@ dataset = xr.open_dataset(demo_lut)

py_wake/tests/test_flow_map.py

@@ -499,1 +499,55 @@ from py_wake.flow_map import HorizontalGrid, YZGrid, Points, XYGrid, XZGrid
     sim_res.flow_map(HorizontalGrid(x=np.linspace(0, 1000, 10), y=np.linspace(0, 1000, 10)), time=[0, 1, 2])
+
+
+def test_wd_dependent_dst():
+    wfm = IEA37CaseStudy1(16)
+    x, y = wfm.site.initial_position.T
+    sim_res = wfm(x, y, wd=np.arange(360), ws=np.arange(3, 25))
+
+    X, Y = np.meshgrid(np.linspace(-1500, 3500, 10), np.linspace(-1500, 1500, 10))
+    dw, hcw = X.flatten(), Y.flatten()
+    dh = dw * 0
+    wf_h = wfm.windTurbines.hub_height()
+    wf_x = 0
+    wf_y = 0
+    from tqdm import tqdm
+    from numpy import newaxis as na
+
+    def run_loop():
+        wd_lst = sim_res.wd.values
+        theta = np.deg2rad(270 - wd_lst)
+        co, si = np.cos(theta), np.sin(theta)
+        x_jl = co[na] * dw[:, na] - hcw[:, na] * si[na] + wf_x
+        y_jl = si[na] * dw[:, na] + hcw[:, na] * co[na] + wf_y
+        h_j = dh + wf_h
+        for wd, x_j, y_j in zip(wd_lst, x_jl.T, y_jl.T):
+            lw_j, WS_eff_jlk, TI_eff_jlk = wfm._flow_map(x_j[:, na], y_j[:, na], h_j[:, na], sim_res.localWind,
+                                                         wd, sim_res.ws, sim_res)
+            if 0:
+                plt.contourf(x_j.reshape(X.shape), y_j.reshape(X.shape),
+                             WS_eff_jlk[:, :, 7].reshape(X.shape), levels=50)
+                wfm.windTurbines.plot(x, y)
+                plt.axis('scaled')
+                plt.show()
+
+    timeit(run_loop, verbose=1, line_profile=0)()
+
+    def run_vec():
+        wd_lst = sim_res.wd.values
+        theta = np.deg2rad(270 - wd_lst)
+        co, si = np.cos(theta), np.sin(theta)
+        x_jl = co[na] * dw[:, na] - hcw[:, na] * si[na] + wf_x
+        y_jl = si[na] * dw[:, na] + hcw[:, na] * co[na] + wf_y
+        h_jl = dh[:, na] + wf_h
+
+        lw_j, WS_eff_jlk, TI_eff_jlk = wfm._flow_map(x_jl, y_jl, h_jl, sim_res.localWind,
+                                                     wd_lst, sim_res.ws, sim_res)
+        if 0:
+            for l, wd_lst in enumerate(wd_lst):
+                plt.contourf(x_jl[:, l].reshape(X.shape), y_jl[:, l].reshape(X.shape),
+                             WS_eff_jlk[:, l, 7].reshape(X.shape), levels=50)
+                wfm.windTurbines.plot(x, y)
+                plt.axis('scaled')
+                plt.show()
+
+    timeit(run_vec, verbose=1, line_profile=0)()

py_wake/tests/test_main.py

 import importlib
 import os
 import pkgutil
+import sys
 import warnings
+from unittest import mock
 
 import pytest
 
-import sys
 import py_wake
-from unittest import mock
 from py_wake.flow_map import Grid
@@ -22,4 +22,7 @@
                 m = importlib.import_module(modname)
+            except DeprecationWarning:
+                pass
             except Exception:
                 print(f"!!!!!! failed to import {modname}")
+
                 raise
@@ -26,0 +29,0 @@

py_wake/tests/test_notebooks.py

@@ -44,3 +44,3 @@ import os
         notebook.remove_empty_end_cell()
-        notebook.check_pip_header()
+        # notebook.check_pip_header()
     except Exception as e:
@@ -47,0 +47,0 @@ raise Exception(notebook.filename + " failed") from e

py_wake/tests/test_sites/test_distances.py

@@ -1,23 +0,31 @@
-from py_wake import np
+import warnings
+
+import matplotlib.pyplot as plt
+import pytest
 from numpy import newaxis as na
-from py_wake.tests import npt
-from py_wake.site.distance import StraightDistance, TerrainFollowingDistance
-from py_wake.site._site import UniformSite
-import pytest
+
+from py_wake import NOJ, np
+from py_wake.deficit_models.gaussian import (
+    BastankhahGaussian,
+    BastankhahGaussianDeficit,
+)
+from py_wake.deficit_models.utils import ct2a_mom1d
+from py_wake.examples.data.hornsrev1 import V80
 from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines
-from py_wake import NOJ
 from py_wake.examples.data.ParqueFicticio import ParqueFicticioSite
-from py_wake.flow_map import HorizontalGrid, XYGrid, XZGrid, Points
-import matplotlib.pyplot as plt
+from py_wake.flow_map import HorizontalGrid, Points, XYGrid, XZGrid
+from py_wake.site._site import UniformSite
+from py_wake.site.distance import StraightDistance, TerrainFollowingDistance
+from py_wake.site.streamline_distance import StreamlineDistance
+from py_wake.tests import npt
+from py_wake.tests.test_wind_farm_models.test_enginering_wind_farm_model import (
+    OperatableV80,
+)
 from py_wake.utils.streamline import VectorField3D
-from py_wake.site.jit_streamline_distance import JITStreamlineDistance
-from py_wake.examples.data.hornsrev1 import V80
+from py_wake.wind_farm_models.engineering_models import (
+    All2AllIterative,
+    PropagateDownwind,
+)
 
-from py_wake.tests.test_wind_farm_models.test_enginering_wind_farm_model import OperatableV80
-from py_wake.wind_farm_models.engineering_models import PropagateDownwind, All2AllIterative
-from py_wake.deficit_models.gaussian import BastankhahGaussianDeficit, BastankhahGaussian
-from py_wake.deficit_models.utils import ct2a_mom1d
-import warnings
 
-
 class FlatSite(UniformSite):
@@ -54,13 +62,16 @@ def __init__(self, distance):
 def test_flat_distances(distance):
-    x = [0, 50, 100, 100, 0]
-    y = [100, 100, 100, 0, 0]
-    h = [0, 10, 20, 30, 0]
-    z = [0, 0, 0, 0]
+    x = np.array([0, 50, 100, 100, 0])
+    y = np.array([100, 100, 100, 0, 0])
+    h = np.array([0, 10, 20, 30, 0])
     wdirs = [-1e-10, 30, 90 + 1e-10]
 
     site = FlatSite(distance=distance)
-    site.distance.setup(src_x_ilk=x, src_y_ilk=y, src_h_ilk=h, src_z_ilk=z)
-    dw_ijlk, hcw_ijlk, dh_ijlk = site.distance(wd_l=np.array(wdirs), WD_ilk=None, src_idx=[0, 1, 2, 3], dst_idx=[4])
-    dw_indices_lkd = site.distance.dw_order_indices(np.array(wdirs))
+    src_idx = [0, 1, 2, 3]
+    dw_ijlk, hcw_ijlk, dh_ijlk = site.distance(src_x_ilk=x[src_idx], src_y_ilk=y[src_idx],
+                                               src_h_ilk=h[src_idx],
+                                               wd_l=np.array(wdirs),
+                                               dst_xyh_jlk=(x[[4]], y[[4]], h[[4]]))
 
+    dw_indices_lkd = site.distance.dw_order_indices(x, y, np.array(wdirs))
+
     if 0:
@@ -78,6 +89,6 @@ distance.plot(wd_ilk=np.array(wdirs)[na, :, na], src_i=[0, 1, 2, 3], dst_i=[4])
                                                      [[-100, -86.6025404, 0]]])
-    npt.assert_array_almost_equal(dh_ijlk[..., 0], [[[0, 0, 0]],
-                                                    [[-10, -10, -10]],
-                                                    [[-20, -20, -20]],
-                                                    [[-30, -30, -30]]])
+    npt.assert_array_almost_equal(dh_ijlk[..., :1, 0], [[[0]],
+                                                        [[-10]],
+                                                        [[-20]],
+                                                        [[-30]]])
     npt.assert_array_equal(dw_indices_lkd[:, 0, :], [[0, 1, 2, 4, 3],
@@ -98,5 +109,6 @@ [2, 1, 0, 3, 4],
     site = FlatSite(distance=distance)
-    site.distance.setup(src_x_ilk=x, src_y_ilk=y, src_h_ilk=h, src_z_ilk=z, dst_xyhz_j=(x, y, [1, 2, 3], z))
-    dw_ijlk, hcw_ijlk, dh_ijlk = site.distance(wd_l=np.array(wdirs), WD_ilk=None)
-    dw_indices_lkd = distance.dw_order_indices(wdirs)
+
+    dw_ijlk, hcw_ijlk, dh_ijlk = site.distance(src_x_ilk=x, src_y_ilk=y, src_h_ilk=h,
+                                               dst_xyh_jlk=(x, y, [1, 2, 3]), wd_l=np.array(wdirs))
+    dw_indices_lkd = distance.dw_order_indices(x, y, wdirs)
     if 0:
@@ -126,3 +138,3 @@ distance.plot(wd_ilk=np.array(wdirs)[na, :, na])
     # check dw indices
-    npt.assert_array_equal(dw_indices_lkd[:, 0], [[1, 0, 2],
+    npt.assert_array_equal(dw_indices_lkd[:, 0], [[0, 1, 2],
                                                   [1, 0, 2]])
@@ -140,4 +152,3 @@
                          ws=site.default_ws)
-    site.distance.setup(x, y, np.zeros_like(x), np.zeros_like(x))
-    dw_iilk, hcw_iilk, _ = site.wt2wt_distances(WD_ilk=lw.WD_ilk, wd_l=None)
+    dw_iilk, hcw_iilk, _ = site.distance(x, y, np.zeros_like(x), WD_ilk=lw.WD_ilk)
     # Wind direction.
@@ -176,3 +187,6 @@ wdir = np.rad2deg(np.arctan2(hcw_iilk, dw_iilk))
     WD = np.array(np.r_[turning, [0] * len(ghost_y)]) + 270
-    sim_res = wfm(np.r_[[0, 500], ghost_x], np.r_[[0, 0], ghost_y], operating=operation, wd=[0, 270], WD=WD)
+    sim_res = wfm(np.r_[[0, 500], ghost_x], np.r_[[0, 0], ghost_y], operating=operation, wd=[270], WD=WD)
+    if 0:
+        sim_res.flow_map(wd=270).plot_wake_map()
+        plt.show()
     with warnings.catch_warnings():
@@ -215,5 +229,5 @@ warnings.simplefilter('ignore', UserWarning)
 
-    hc.distance.setup(src_x_ilk=src_x, src_y_ilk=src_y, src_h_ilk=src_x * 0, src_z_ilk=src_x * 0,
-                      dst_xyhz_j=(dst_x, dst_y, dst_x * 0, dst_x * 0))
-    dw_ijlk, hcw_ijlk, _ = hc.distance(wd_l=np.array([0, 90]), WD_ilk=None)
+    dw_ijlk, hcw_ijlk, _ = hc.distance(src_x_ilk=src_x, src_y_ilk=src_y, src_h_ilk=src_x * 0,
+                                       wd_l=np.array([0, 90]), WD_ilk=None,
+                                       dst_xyh_jlk=(dst_x, dst_y, dst_x * 0))
 
@@ -273,4 +287,3 @@ if 0:
     y_j = x_j * 0
-    site.distance.setup([-200], [0], [130], [0], (x_j, y_j, y_j + 130, y_j * 0))
-    d = site.distance(wd_l=[270])[0][0, :, 0, 0]
+    d = site.distance([-200], [0], [130], [0], wd_l=[270], dst_xyh_jlk=(x_j, y_j, y_j + 130))[0][0, :, 0, 0]
 
@@ -294,10 +307,12 @@ ref = x_j - x_j[0]
 
-    x = [0, 50, 100, 100, 0]
-    y = [100, 100, 100, 0, 0]
-    h = [0, 10, 20, 30, 0]
-    z = [0, 0, 0, 0, 0]
+    x = np.array([0, 50, 100, 100, 0])
+    y = np.array([100, 100, 100, 0, 0])
+    h = np.array([0, 10, 20, 30, 0])
     wdirs = [0, 30, 90]
     distance = StraightDistance()
-    distance.setup(src_x_ilk=x, src_y_ilk=y, src_h_ilk=h, src_z_ilk=z)
-    distance.plot(wd_l=np.array(wdirs), src_idx=[0], dst_idx=[3])
+
+    src_idx = [0]
+    dst_idx = [3]
+    distance.plot(src_x_ilk=x[src_idx], src_y_ilk=y[src_idx], src_h_ilk=h[src_idx], wd_l=np.array(wdirs),
+                  dst_xyh_jlk=[x[dst_idx], y[dst_idx], h[dst_idx]])
     if 0:
@@ -308,7 +323,7 @@ plt.show()
 
-def test_JITStreamlinesparquefictio():
+def test_Streamlinesparquefictio():
     site = ParqueFicticioSite()
     wt = IEA37_WindTurbines()
     vf3d = VectorField3D.from_WaspGridSite(site)
-    site.distance = JITStreamlineDistance(vf3d)
+    site.distance = StreamlineDistance(vf3d)
 
@@ -319,3 +334,3 @@ x, y = site.initial_position[:].T
     sim_res = wfm(x, y, wd=wd, ws=10)
-    dw = site.distance(wd_l=wd, WD_ilk=np.repeat(wd[na, na], len(x), 0))[0][:, :, 0, 0]
+    dw = site.distance(x, y, x * 0 + wt.hub_height(), wd_l=wd, WD_ilk=np.repeat(wd[na, na], len(x), 0))[0][:, :, 0, 0]
     # streamline downwind distance (positive numbers, upper triangle) cannot be shorter than
@@ -338,3 +353,3 @@ # straight line distances in opposite direction (negative numbers, lower triangle)
 
-def test_JITStreamlinesparquefictio_yz():
+def test_Streamlinesparquefictio_yz():
     site = ParqueFicticioSite()
@@ -345,3 +360,3 @@ site.ds.Turning[:] *= 0
     vf3d = VectorField3D.from_WaspGridSite(site)
-    site.distance = JITStreamlineDistance(vf3d)
+    site.distance = StreamlineDistance(vf3d)
 
@@ -356,3 +371,4 @@ x, y = site.initial_position[3].T
     sim_res = wfm(wt_x, wt_y, wd=wd, ws=10)
-    dw = site.distance(wd_l=wd, WD_ilk=np.repeat(wd[na, na], len(wt_x), 0))[0][:, :, 0, 0]
+    dw = site.distance(wt_x, wt_y, wt_x * 0 + wt.hub_height(), wd_l=wd,
+                       WD_ilk=np.repeat(wd[na, na], len(wt_x), 0))[0][:, :, 0, 0]
     # streamline downwind distance (positive numbers, upper triangle) cannot be shorter than
@@ -359,0 +375,0 @@ # straight line distances in opposite direction (negative numbers, lower triangle)

py_wake/tests/test_sites/test_site.py

@@ -220,4 +220,3 @@ import os
                          ws=site.default_ws)
-    site.distance.setup(x, y, np.zeros_like(x), np.zeros_like(x))
-    dw_iil, hcw_iil, _ = site.wt2wt_distances(wd_l=lw.wd)
+    dw_iil, hcw_iil, _ = site.distance(x, y, x * 0, wd_l=lw.wd)
     # Wind direction.
@@ -227,3 +226,3 @@ wdir = np.rad2deg(np.arctan2(hcw_iil, dw_iil))
         wdir[:, 0, 0, 0],
-        [180, -90, -18, 54, 126, -162, -90, -54, -18, 18, 54, 90, 126, 162, -162, -126],
+        [0, -90, -18, 54, 126, -162, -90, -54, -18, 18, 54, 90, 126, 162, -162, -126],
         atol=1e-4)
@@ -230,0 +229,0 @@

py_wake/tests/test_sites/test_wasp_grid_site.py

@@ -203,5 +203,5 @@ from py_wake import np
     x, y = site.initial_position.T
-    site.distance.setup(src_x_ilk=x, src_y_ilk=y, src_h_ilk=np.array([70]), src_z_ilk=x * 0,
-                        dst_xyhz_j=(x, y, np.array([70]), x * 0))
-    dw_ijlk, cw_ijlk, dh_ijlk = site.distance(wd_l=[0])
+    dw_ijlk, cw_ijlk, dh_ijlk = site.distance(src_x_ilk=x, src_y_ilk=y, src_h_ilk=np.array([70]), wd_l=[0],
+                                              dst_xyh_jlk=(x, y, np.array([70])))
+
     npt.assert_almost_equal(dw_ijlk[0, :, 0, 0], dw_ref)
@@ -208,0 +208,0 @@

py_wake/tests/test_utils/test_gradient_utils.py

@@ -11,3 +11,4 @@ from autograd import numpy as anp
 from py_wake.utils import gradients
-from py_wake.utils.gradients import autograd, plot_gradients, fd, cs, hypot, cabs, interp, set_gradient_function
+from py_wake.utils.gradients import autograd, plot_gradients, fd, cs, hypot, cabs, interp, set_gradient_function, \
+    item_assign
 from py_wake.wind_turbines import WindTurbines
@@ -509,1 +510,22 @@ from py_wake.wind_turbines import _wind_turbines
     autograd(f)([2, 3, 4])
+
+
[e-mail]('idx', [[0, 3, 4], slice(1, 6, 2)])
[e-mail]('axis', [0, 1])
+def test_item_assign(idx, axis):
+
+    def f(x):
+        if axis == 0:
+            v = x[idx]
+        else:
+            v = x[:, idx]
+        x = item_assign(x + 0., idx=idx, values=v * 2, axis=axis)
+        return x.sum()
+
+    x = np.arange(56).reshape((8, 7))
+    v = x.take(np.arange(x.shape[axis])[idx], axis) * 2
+    npt.assert_array_equal(item_assign(x + 0., idx, values=v, axis=axis),
+                           item_assign(x + 0., idx, values=v, axis=axis, test=1))
+
+    npt.assert_array_almost_equal(fd(f)(x), cs(f)(x))
+    npt.assert_array_equal(cs(f)(x), autograd(f)(x))

py_wake/tests/test_utils/test_model_gradients.py

@@ -1,32 +0,47 @@
-from numpy import newaxis as na
+import warnings
+
+import matplotlib.pyplot as plt
 import pytest
 import xarray as xr
+from numpy import newaxis as na
 
-import matplotlib.pyplot as plt
 from py_wake import np
-from py_wake.deficit_models.deficit_model import WakeDeficitModel, BlockageDeficitModel
-from py_wake.deficit_models.gaussian import BastankhahGaussianDeficit, NiayifarGaussianDeficit
+from py_wake.deficit_models.deficit_model import BlockageDeficitModel, WakeDeficitModel
+from py_wake.deficit_models.gaussian import (
+    BastankhahGaussianDeficit,
+    NiayifarGaussianDeficit,
+)
 from py_wake.deficit_models.no_wake import NoWakeDeficit
+from py_wake.deficit_models.noj import NOJDeficit
 from py_wake.deflection_models.deflection_model import DeflectionModel
-from py_wake.examples.data.ParqueFicticio._parque_ficticio import ParqueFicticioSite
 from py_wake.examples.data.hornsrev1 import V80, Hornsrev1Site
 from py_wake.examples.data.iea34_130rwt._iea34_130rwt import IEA34_130_1WT_Surrogate
-from py_wake.examples.data.iea37._iea37 import IEA37_WindTurbines, IEA37Site, IEA37WindTurbines
+from py_wake.examples.data.iea37._iea37 import (
+    IEA37_WindTurbines,
+    IEA37Site,
+    IEA37WindTurbines,
+)
+from py_wake.examples.data.ParqueFicticio._parque_ficticio import ParqueFicticioSite
 from py_wake.ground_models.ground_models import GroundModel
+from py_wake.rotor_avg_models.area_overlap_model import AreaOverlapAvgModel
 from py_wake.rotor_avg_models.rotor_avg_model import RotorAvgModel
 from py_wake.site.distance import StraightDistance
-from py_wake.site.shear import LogShear, PowerShear, Shear, MOSTShear
-from py_wake.superposition_models import SuperpositionModel, AddedTurbulenceSuperpositionModel
+from py_wake.site.shear import LogShear, MOSTShear, PowerShear, Shear
+from py_wake.site.streamline_distance import StreamlineDistance
+from py_wake.site.xrsite import XRSite
+from py_wake.superposition_models import (
+    AddedTurbulenceSuperpositionModel,
+    SuperpositionModel,
+)
 from py_wake.tests import npt
-from py_wake.turbulence_models.stf import STF2017TurbulenceModel, STF2005TurbulenceModel
+from py_wake.turbulence_models.stf import STF2005TurbulenceModel, STF2017TurbulenceModel
 from py_wake.turbulence_models.turbulence_model import TurbulenceModel
 from py_wake.utils import gradients
-from py_wake.utils.gradients import autograd, plot_gradients, fd, cs
+from py_wake.utils.gradients import autograd, cs, fd, plot_gradients
 from py_wake.utils.model_utils import get_models
-from py_wake.wind_farm_models.engineering_models import PropagateDownwind, All2AllIterative, EngineeringWindFarmModel
-from py_wake.deficit_models.noj import NOJDeficit
-from py_wake.site.jit_streamline_distance import JITStreamlineDistance
-from py_wake.site.xrsite import XRSite
-from py_wake.rotor_avg_models.area_overlap_model import AreaOverlapAvgModel
-import warnings
+from py_wake.wind_farm_models.engineering_models import (
+    All2AllIterative,
+    EngineeringWindFarmModel,
+    PropagateDownwind,
+)
 
@@ -239,3 +254,3 @@
 def test_distance_models(model):
-    if model not in [None, JITStreamlineDistance]:
+    if model not in [None, StreamlineDistance]:
         site = ParqueFicticioSite(distance=model())
@@ -242,0 +257,0 @@ x, y = site.initial_position[3]

py_wake/tests/test_utils/test_profiling.py

@@ -20,6 +20,6 @@ import time
         try:
-            # check 3 successive runs (all should take .1 s and 10MB
+            # check 3 successive runs (all should take 1 s and 10MB
             for _ in range(3):
-                r, t, m = profileit(f)(.1, 100)
-                npt.assert_allclose(t, 0.1, rtol=.2)
+                r, t, m = profileit(f)(1, 100)
+                npt.assert_allclose(t, 1, rtol=.2)
                 npt.assert_allclose(m, 100, rtol=.2)
@@ -26,0 +26,0 @@ return

py_wake/tests/test_wind_farm_models/test_parallel_daep.py

@@ -46,2 +46,3 @@ import numpy.testing as npt
     wfm = Nygaard_2022(site, V80())
+    wd = ds.WD.mean(['x', 'y', 'h'])
 
@@ -60,3 +61,3 @@ test_x, test_y = [
         _f = np.zeros_like(T)
-        sim_res = wfm(x, y, time=T, ws=_f, wd=_f, n_cpu=n_cpu)
+        sim_res = wfm(x, y, time=T, ws=_f, wd=wd, n_cpu=n_cpu)
         return sim_res.aep().sum().values
@@ -72,3 +73,3 @@
             ws=_f,
-            wd=_f,
+            wd=wd,
             time=T,
@@ -75,0 +76,0 @@ n_cpu=n_cpu,

py_wake/utils/gradients.py

@@ -39,2 +39,28 @@ from py_wake import np
 
+def item_assign(x, idx, values, axis=0, test=0):
+    if axis:
+        x = np.moveaxis(x, axis, 0)
+        values = np.moveaxis(values, axis, 0)
+    assert x[idx].shape == np.shape(values)
+    if isinstance(x, ArrayBox) or test:
+        if isinstance(idx, slice):
+            idx = np.arange(x.shape[axis])[idx]
+        res = []
+        idx_i = 0
+        for x_i in range(x.shape[0]):
+            if idx_i < len(idx) and x_i == idx[idx_i]:
+                res.append(values[idx_i])
+                idx_i += 1
+            else:
+                res.append(x[x_i])
+        x = np.array(res)
+
+    else:
+        # x = x.copy()
+        x[idx] = values
+    if axis:
+        x = np.moveaxis(x, 0, axis)
+    return x
+
+
 # def asanyarray(x, dtype=None, order=None):
@@ -41,0 +67,0 @@ #     if isinstance(x, (ArrayBox)):

py_wake/utils/model_utils.py

@@ -176,3 +176,3 @@ import inspect
     from py_wake.ground_models.ground_models import NoGround
-    from py_wake.site.jit_streamline_distance import JITStreamlineDistance
+    from py_wake.site.streamline_distance import StreamlineDistance
     return {
@@ -193,3 +193,3 @@ "WindFarmModel": ([EngineeringWindFarmModel], [], PropagateDownwind),
         "Shear": ([], [], None),
-        "StraightDistance": ([], [JITStreamlineDistance], None),
+        "StraightDistance": ([], [StreamlineDistance], None),
 
@@ -287,5 +287,5 @@ }
     x, y = map(np.asarray, [x, y])
-    wfm.site.distance.setup(src_x_ilk=[[[0]]], src_y_ilk=[[[0]]], src_h_ilk=[[[0]]], src_z_ilk=[[[0]]],
-                            dst_xyhz_j=(x, y, x * 0, x * 0))
-    dw_ijlk, hcw_ijlk, dh_ijlk = wfm.site.distance(wd_l=wd)
+
+    dw_ijlk, hcw_ijlk, dh_ijlk = wfm.site.distance(src_x_ilk=[[[0]]], src_y_ilk=[[[0]]], src_h_ilk=[[[0]]],
+                                                   dst_xyh_jlk=(x, y, x * 0), wd_l=wd)
     sim_res = wfm([0], [0], ws=ws, wd=wd, **kwargs)
@@ -292,0 +292,0 @@

py_wake/utils/parallelization.py

@@ -55,14 +55,14 @@ import multiprocessing
 
-def get_map_func(n_cpu, verbose, desc='', unit='it'):
-    n_cpu = n_cpu or multiprocessing.cpu_count()
-    if n_cpu > 1:
-        map_func = get_pool_map(n_cpu)
-    else:
-        from tqdm import tqdm
+# def get_map_func(n_cpu, verbose, desc='', unit='it'):
+#     n_cpu = n_cpu or multiprocessing.cpu_count()
+#     if n_cpu > 1:
+#         map_func = get_pool_map(n_cpu)
+#     else:
+#         from tqdm import tqdm
+#
+#         def map_func(f, iter):
+#             return tqdm(map(f, iter), desc=desc, unit=unit, total=len(iter), disable=not verbose)
+#     return map_func
 
-        def map_func(f, iter):
-            return tqdm(map(f, iter), desc=desc, unit=unit, total=len(iter), disable=not verbose)
-    return map_func
 
-
 def get_starmap_func(n_cpu, verbose, desc='', unit='it', leave=True):
@@ -69,0 +69,0 @@ n_cpu = n_cpu or multiprocessing.cpu_count()

py_wake/utils/rans_lut_utils.py

@@ -216,3 +216,4 @@ import xarray as xr
             setattr(sim_res_wd, k, getattr(sim_res, k))
-        WS_eff_jlk = np.concatenate([sim_res.windFarmModel._flow_map(*pos_j.T, sim_res_wd)[1]
+        wd, ws = sim_res_wd.wd.values, sim_res_wd.ws.values
+        WS_eff_jlk = np.concatenate([sim_res.windFarmModel._flow_map(*[p[:, na] for p in pos_j.T], sim_res.localWind, wd, ws, sim_res_wd)[1]
                                      for pos_j in np.array_split(cxyzg_xj.T, nflowmaps)])
@@ -224,2 +225,20 @@
 
+    # Vectorized version of get_WS_eff_ik.
+    # Against expectations, this version is slower with numpy 2.3.4 and scipy 1.16.3.
+    # The slowdown seems to
+    # def get_WS_eff_ilk():
+    #     # project rotorAvgModel node positions to plane perpendicular to current wind direction
+    #     o3 = o3_l[:, :]
+    #     o1 = np.cross(o2, o3.T).T
+    #     cxyzg_xijl = nx_in[na, :, :, na] * o3[:, na, na, :, ] + ny_in[na, :, :, na] * \
+    #         o1[:, na, na, :] + nz_in[na, :, :, na] * o2[:, na, na, na] + pos[:, :, na, na]
+    #     cxyzg_xjl = cxyzg_xijl.reshape((3, -1, L))
+    #     wd, ws = sim_res.wd.values, sim_res.ws.values
+    #     WS_eff_jlk = sim_res.windFarmModel._flow_map(*cxyzg_xjl, sim_res.localWind, wd, ws, sim_res)[1]
+    #
+    #     # Integrate U over AD
+    #     return np.sum(np.reshape(WS_eff_jlk, (I, ntheta * (nr - 1), L, K)) * weights[:, :, na], axis=1)
+    #
+    # WS_eff_star_ilk = get_WS_eff_ilk()
+
     U_tab_lst = [U_CT_CP_AD[0, :] for U_CT_CP_AD in aDControl.U_CT_CP_AD]
@@ -226,0 +245,0 @@

py_wake/utils/streamline.py

@@ -1,7 +0,8 @@
-from py_wake.utils.grid_interpolator import GridInterpolator
-from py_wake import np
+import xarray as xr
 from numpy import newaxis as na
-import xarray as xr
 
+from py_wake import np
+from py_wake.utils.grid_interpolator import GridInterpolator
 
+
 class VectorField3D():
@@ -20,3 +21,3 @@ def __init__(self, da):
 
-    def __call__(self, wd, x, y, h):
+    def __call__(self, wd, time, x, y, h):
         return self.interpolator(np.array([np.atleast_1d(v) for v in [wd, x, y, h]]).T, bounds='limit')
@@ -34,4 +35,13 @@
 
-    def stream_lines(self, wd, start_points, dw_stop, step_size=20):
-        wd = np.full(len(start_points), wd)
+    def stream_lines(self, wd, start_points, dw_stop, time=None, step_size=20):
+        # print(np.shape(wd), np.shape(start_points), time)
+        if time is None:
+            wd = np.full(len(start_points), wd)
+            time = wd * 0
+        # else:
+        #     if len(start_points) == 1:
+        #         # 1 start point, multiple wd/time
+        #         start_points = np.broadcast_to(start_points, (len(wd), 3))
+        #     elif len(time) == 1:
+        #         time = wd * 0 + time
         stream_lines = [start_points]
@@ -42,3 +52,3 @@ m = np.arange(len(wd))
             p = stream_lines[-1].copy()
-            v = self(wd[m], p[m, 0], p[m, 1], p[m, 2])
+            v = self(wd[m], time[m], p[m, 0], p[m, 1], p[m, 2])
             v = v / (np.sqrt(np.sum(v**2, -1)) / step_size)[:, na]  # normalize vector distance to step_size
@@ -45,0 +55,0 @@ p[m] += v

py_wake/version.py

@@ -31,5 +31,5 @@ # file generated by setuptools-scm
 
-__version__ = version = '2.6.13'
-__version_tuple__ = version_tuple = (2, 6, 13)
+__version__ = version = '2.6.14'
+__version_tuple__ = version_tuple = (2, 6, 14)
 
 __commit_id__ = commit_id = None

py_wake/wind_farm_models/engineering_models.py

@@ -7,3 +7,3 @@ from abc import abstractmethod
 from py_wake.deflection_models.deflection_model import DeflectionModel
-from py_wake.utils.gradients import cabs
+from py_wake.utils.gradients import cabs, item_assign
 from py_wake.rotor_avg_models.rotor_avg_model import RotorAvgModel, NodeRotorAvgModel, WSPowerRotorAvg
@@ -19,4 +19,5 @@ from py_wake.turbulence_models.turbulence_model import TurbulenceModel
 from py_wake.deficit_models.no_wake import NoWakeDeficit
-from py_wake.utils.parallelization import get_map_func, get_starmap_func
+from py_wake.utils.parallelization import get_starmap_func
 import multiprocessing
+from py_wake.wind_farm_models.external_wind_farm_models import ExternalWindFarm
 
@@ -54,3 +55,4 @@
     def __init__(self, site, windTurbines: WindTurbines, wake_deficitModel, superpositionModel, rotorAvgModel=None,
-                 blockage_deficitModel=None, deflectionModel=None, turbulenceModel=None, inputModifierModels=[]):
+                 blockage_deficitModel=None, deflectionModel=None, turbulenceModel=None, inputModifierModels=[],
+                 externalWindFarms=[]):
 
@@ -65,3 +67,4 @@ WindFarmModel.__init__(self, site, windTurbines)
                                   (turbulenceModel, TurbulenceModel, 'turbulenceModel')] +
-                                 [(imm, InputModifierModel, 'inputModificierModels') for imm in inputModifierModels]):
+                                 [(imm, InputModifierModel, 'inputModificierModels') for imm in inputModifierModels] +
+                                 [(ewf, ExternalWindFarm, 'externalWindFarms') for ewf in externalWindFarms]):
             check_model(model, cls, name)
@@ -72,2 +75,3 @@ if model is not None:
         self.inputModifierModels = inputModifierModels
+        self.externalWindFarms = externalWindFarms
 
@@ -81,2 +85,4 @@ if isinstance(superpositionModel, (WeightedSum, CumulativeWakeSum)):
                 'WeightedSum and CumulativeWakeSum does not work WSPowerRotorAvg'
+            assert len(externalWindFarms) == 0, \
+                'WeightedSum and CumulativeWakeSum does not work with ExternalWindFarms'
         # TI_eff requires a turbulence model
@@ -123,2 +129,4 @@ assert 'TI_eff_ilk' not in wake_deficitModel.args4deficit or turbulenceModel
             self.args4all |= set(input_modifier.args4model)
+        for ext_wf in self.externalWindFarms:
+            self.args4all |= set(ext_wf.args4model)
 
@@ -238,2 +246,4 @@ def __str__(self):
             kwargs['TI_eff_ilk'] = lw.TI_ilk + 0.  # autograd-friendly copy
+        if 'time' in lw:
+            kwargs['time'] = lw.time
 
@@ -265,6 +275,2 @@ self._check_input(kwargs)
         # Calculate down-wind and cross-wind distances
-        z_ilk = self.site.elevation(kwargs['x_ilk'], kwargs['y_ilk'])
-        self.site.distance.setup(
-            np.asarray(kwargs['x_ilk']), np.asarray(kwargs['y_ilk']), np.asarray(kwargs['h_ilk']), np.asarray(z_ilk))
-
         WS_eff_ilk, TI_eff_ilk, ct_ilk, kwargs = self._calc_wt_interaction(**kwargs)
@@ -280,10 +286,13 @@
 
-    def get_map_args(self, x_j, y_j, h_j, sim_res_data, D_dst=0):
+    def get_map_args(self, x_j, y_j, h_j, lw, wd, ws, sim_res_data, D_dst=0):
         type_i = sim_res_data.type.values
         wt_d_i = self.windTurbines.diameter(type_i)
-        wd, ws = [np.atleast_1d(sim_res_data[k].values) for k in ['wd', 'ws']]
         if 'time' in sim_res_data:
             time = np.atleast_1d(sim_res_data['time'].values)
+            map_arg_funcs = {'time': lambda l, j: time[l]}
         else:
             time = False
+            map_arg_funcs = {}
+        if 'wd' in sim_res_data.dims:
+            sim_res_data = sim_res_data.sel(wd=wd)
         wt_x_ilk = sim_res_data['x'].ilk()
@@ -300,16 +309,17 @@ WD_il = sim_res_data.WD.ilk()
 
-            def wrap(l):
+            def wrap(l, j):
                 l_ = [l, slice(0, 1)][v.shape[1] == 1]
                 return v[:, l_]
             return wrap
-        map_arg_funcs = {k.replace('CT', 'ct') + '_ilk': get_ilk(k)
-                         for k in sim_res_data if k not in ['wd_bin_size', 'ws_l', 'ws_u']}
+        map_arg_funcs.update({k.replace('CT', 'ct') + '_ilk': get_ilk(k)
+                              for k in sim_res_data if k not in ['wd_bin_size', 'ws_l', 'ws_u']})
         map_arg_funcs.update({
-            'type_il': lambda l: type_i[:, na],
-            'D_src_il': lambda l: wt_d_i[:, na],
-            'D_dst_ijl': lambda l: np.zeros((1, 1, 1)) + D_dst,
-            'IJLK': lambda l=slice(None), I=I, J=J, L=L, K=K: (I, J, len(np.arange(L)[l]), K)})
+            'type_il': lambda l, j: type_i[:, na],
+            'D_src_il': lambda l, j: wt_d_i[:, na],
+            'D_dst_ijl': lambda l, j: np.zeros((1, 1, 1)) + D_dst,
+            'IJLK': lambda l=slice(None), j=slice(None), I=I, J=J, L=L, K=K: (I, len(np.arange(J)[j]), len(np.arange(L)[l]), K)})
+
         for k in ['WS', 'WD', 'TI']:
             if k in sim_res_data:
-                if k + '_ilk' in sim_res_data.localWind.overwritten:
+                if k + '_ilk' in lw.overwritten:
                     if 'wt' not in sim_res_data[k].dims and 'i' not in sim_res_data[k].dims:
@@ -320,9 +330,7 @@ lw_j.add_ilk(k + "_ilk", sim_res_data[k])
                             f"The WT dependent {k} that was provided for the simulation is not available at the flow map points and therefore ignored")
-        return map_arg_funcs, lw_j, wd, WD_il
+        return map_arg_funcs, lw_j, WD_il
 
-    def _get_flow_l(self, model_kwargs, wt_x_ilk, wt_y_ilk, wt_h_ilk, x_j, y_j, h_j, wd, WD_ilk, WS_jlk, TI_jlk):
-        wt_z_ilk = self.site.elevation(wt_x_ilk, wt_y_ilk)
-        z_j = self.site.elevation(x_j, y_j)
-        self.site.distance.setup(wt_x_ilk, wt_y_ilk, wt_h_ilk, wt_z_ilk, (x_j, y_j, h_j, z_j))
-        dw_ijlk, hcw_ijlk, dh_ijlk = self.site.distance(wd_l=wd, WD_ilk=WD_ilk)
+    def _get_flow_l(self, model_kwargs, wt_x_ilk, wt_y_ilk, wt_h_ilk, x_jl, y_jl, h_jl, wd, WD_ilk, WS_jlk, TI_jlk):
+        dw_ijlk, hcw_ijlk, dh_ijlk = self.site.distance(wt_x_ilk, wt_y_ilk, wt_h_ilk, wd_l=wd, WD_ilk=WD_ilk,
+                                                        time=model_kwargs.get('time', None), dst_xyh_jlk=(x_jl, y_jl, h_jl))
 
@@ -362,2 +370,10 @@ if self.wec != 1:
         # ===============================================================================================================
+        # Replace deficit from external farms
+        # ===============================================================================================================
+        for i, ewf in enumerate(self.externalWindFarms, len(wt_x_ilk) - len(self.externalWindFarms)):
+            deficit_ijlk[i] = ewf(i=i, l=np.ones(len(wd), dtype=bool), deficit_jlk=None,
+                                  **model_kwargs,
+                                  dst_xyh_jlk=[x_jl[:, :, na], y_jl[:, :, na], h_jl[:, :, na]]
+                                  )
+        # ===============================================================================================================
         # Calculate added Turbulence
@@ -397,4 +413,5 @@ # ===============================================================================================================
 
-    def _aep_map(self, x_j, y_j, h_j, type_j, sim_res_data, memory_GB=1, n_cpu=1):
-        lw_j, WS_eff_jlk, _ = self._flow_map(x_j, y_j, h_j, sim_res_data, memory_GB=memory_GB, n_cpu=n_cpu)
+    def _aep_map(self, x_jl, y_jl, h_jl, type_j, sim_res_data, memory_GB=1, n_cpu=1):
+        lw_j, WS_eff_jlk, _ = self._flow_map(x_jl, y_jl, h_jl, sim_res_data.localWind, sim_res_data.wd.values,
+                                             sim_res_data.ws.values, sim_res_data, memory_GB=memory_GB, n_cpu=n_cpu)
         power_kwargs = {}
@@ -409,10 +426,11 @@ if 'type' in (self.windTurbines.powerCtFunction.required_inputs +
 
-    def _flow_map(self, x_j, y_j, h_j, sim_res_data, D_dst=0, memory_GB=1, n_cpu=1):
+    def _flow_map(self, x_jl, y_jl, h_jl, lw, wd, ws, sim_res_data, D_dst=0, memory_GB=1, n_cpu=1):
         """call this function via SimulationResult.flow_map"""
-        arg_funcs, lw_j, wd, WD_il = self.get_map_args(x_j, y_j, h_j, sim_res_data, D_dst=D_dst)
+        wd, ws = np.atleast_1d(wd), np.atleast_1d(ws)
+        arg_funcs, lw_j, WD_il = self.get_map_args(x_jl, y_jl, h_jl, lw, wd, ws, sim_res_data, D_dst=D_dst)
         I, J, L, K = arg_funcs['IJLK']()
 
         if I == 0:
-            return (lw_j, np.broadcast_to(lw_j.WS_ilk, (len(x_j), L, K)).astype(float),
-                    np.broadcast_to(lw_j.TI_ilk, (len(x_j), L, K)).astype(float))
+            return (lw_j, np.broadcast_to(lw_j.WS_ilk, (len(x_jl), L, K)).astype(float),
+                    np.broadcast_to(lw_j.TI_ilk, (len(x_jl), L, K)).astype(float))
         n_cpu = n_cpu or multiprocessing.cpu_count()
@@ -430,13 +448,16 @@ # *6=dx_ijlk, dy_ijlk, dz_ijlk, dh_ijlk, deficit, blockage
 
-        xyh_j_lst = list(zip(*[np.array_split(x_j, j_chunks),
-                               np.array_split(y_j, j_chunks),
-                               np.array_split(h_j, j_chunks)]))
+        def get_jl_args(j, l):
+            js, ls = slice(j[0], j[-1] + 1), slice(l[0], l[-1] + 1)
+            if x_jl.shape[1] == 1:
+                ls = slice(None)
 
-        l_iter = [({k: arg_funcs[k](l) for k in arg_funcs},
-                   *[(v, v[:, slice(l[0], l[-1] + 1)])[np.shape(v)[1] == L] for v in [wt_x_ilk, wt_y_ilk, wt_h_ilk]],
-                   x_j_, y_j_, h_j_, wd[l], WD_il[:, l],
-                   lw_j.WS_ilk[:, (l, [0])[lw_j.WS_ilk.shape[1] == 1]],
-                   lw_j.TI_ilk[:, (l, [0])[lw_j.TI_ilk.shape[1] == 1]])
+            return ({k: arg_funcs[k](l, j) for k in arg_funcs},
+                    *[(v, v[:, slice(l[0], l[-1] + 1)])[np.shape(v)[1] == L] for v in [wt_x_ilk, wt_y_ilk, wt_h_ilk]],
+                    x_jl[js, ls], y_jl[js, ls], np.broadcast_to(h_jl, x_jl.shape)[js, ls], wd[l], WD_il[:, l],
+                    lw_j.WS_ilk[:, (l, [0])[lw_j.WS_ilk.shape[1] == 1]],
+                    lw_j.TI_ilk[:, (l, [0])[lw_j.TI_ilk.shape[1] == 1]])
+
+        l_iter = [get_jl_args(j, l)
                   for l in np.array_split(np.arange(L).astype(int), wd_chunks)
-                  for x_j_, y_j_, h_j_ in xyh_j_lst]
+                  for j in np.array_split(np.arange(J).astype(int), j_chunks)]
         WS_eff_jlk, TI_eff_jlk = zip(*map_func(self._get_flow_l, l_iter))
@@ -526,3 +547,3 @@ WS_eff_jlk = np.concatenate([np.concatenate(WS_eff_jlk[i * j_chunks:(i + 1) * j_chunks], 0)
         I = kwargs['I']
-        dw_order_indices_ld = self.site.distance.dw_order_indices(wd)[:, 0]
+        dw_order_indices_ld = self.site.distance.dw_order_indices(kwargs['x_ilk'], kwargs['y_ilk'], wd)[:, 0]
 
@@ -732,4 +753,8 @@ if self.blockage_deficitModel:
 
+                def get_pos(v_ilk, i_wt_l, i_wd_l):
+                    return v_ilk[i_wt_l, np.minimum(i_wd_l, v_ilk.shape[1] - 1).astype(int)]
+                dst_xyh_jlk = [get_pos(kwargs[k + '_ilk'], i_dw.T, i_wd_l) for k in 'xyh']
                 dw_ijlk, hcw_ijlk, dh_ijlk = self.site.distance(
-                    wd_l=wd, WD_ilk=WD_mk[m][na], src_idx=i_wt_l, dst_idx=i_dw.T)
+                    *[get_pos(kwargs[k + '_ilk'], i_wt_l, i_wd_l)[na] for k in 'xyh'], wd_l=wd, WD_ilk=WD_mk[m][na],
+                    dst_xyh_jlk=dst_xyh_jlk)
 
@@ -787,2 +812,13 @@ for inputModidifierModel in self.inputModifierModels:
                     deficit, blockage = self._calc_deficit(**model_kwargs)
+
+                    for i, ewf in enumerate(self.externalWindFarms, I - len(self.externalWindFarms)):
+                        # print(i)
+                        deficit = ewf(i=0, l=i_wt_l == i,
+                                      deficit_jlk=deficit[0],
+
+                                      **model_kwargs,
+                                      # **{**model_kwargs, 'dw_ijlk': sdw_ijlk, 'hcw_ijlk': shcw_ijlk, 'dh_ijlk': sdh_ijlk}
+                                      dst_xyh_jlk=dst_xyh_jlk
+                                      )[na]
+
                     if self.blockage_deficitModel:
@@ -822,3 +858,4 @@ blockage_nk.append(blockage[0])
                  deflectionModel=None, turbulenceModel=None, rotorAvgModel=None,
-                 inputModifierModels=[]):
+                 inputModifierModels=[],
+                 externalWindFarms=[]):
         """Initialize flow model
@@ -847,3 +884,4 @@
                                           blockage_deficitModel=None, deflectionModel=deflectionModel,
-                                          turbulenceModel=turbulenceModel, inputModifierModels=inputModifierModels)
+                                          turbulenceModel=turbulenceModel, inputModifierModels=inputModifierModels,
+                                          externalWindFarms=externalWindFarms)
         self.direction = 'down'
@@ -857,3 +895,3 @@
 
-        dw_order_indices_ld = self.site.distance.dw_order_indices(wd)[:, 0]
+        dw_order_indices_ld = self.site.distance.dw_order_indices(kwargs['x_ilk'], kwargs['y_ilk'], wd)[:, 0]
         return self._propagate_deficit(wd, dw_order_indices_ld, WS_ilk, **kwargs)
@@ -869,3 +907,3 @@
                  blockage_deficitModel=None, deflectionModel=None, turbulenceModel=None,
-                 convergence_tolerance=1e-6, rotorAvgModel=None, inputModifierModels=[]):
+                 convergence_tolerance=1e-6, rotorAvgModel=None, inputModifierModels=[], externalWindFarms=[]):
         """Initialize flow model
@@ -901,3 +939,4 @@
                                           blockage_deficitModel=blockage_deficitModel, deflectionModel=deflectionModel,
-                                          turbulenceModel=turbulenceModel, inputModifierModels=inputModifierModels)
+                                          turbulenceModel=turbulenceModel, inputModifierModels=inputModifierModels,
+                                          externalWindFarms=externalWindFarms)
         self.convergence_tolerance = convergence_tolerance
@@ -921,3 +960,3 @@
             self.blockage_deficitModel = None
-            dw_order_indices_ld = self.site.distance.dw_order_indices(wd)[:, 0]
+            dw_order_indices_ld = self.site.distance.dw_order_indices(kwargs['x_ilk'], kwargs['y_ilk'], wd)[:, 0]
             self.direction = 'down'
@@ -939,3 +978,5 @@ WS_eff_ilk = PropagateUpDownIterative._propagate_deficit(
         a = 1
-        dw_iilk, hcw_iilk, dh_iilk = self.site.distance(wd_l=wd, WD_ilk=WD_ilk)
+
+        dst_xyh_jlk = [kwargs[k + '_ilk'] for k in 'xyh']
+        dw_iilk, hcw_iilk, dh_iilk = self.site.distance(*[kwargs[k + '_ilk'] for k in 'xyh'], wd_l=wd, WD_ilk=WD_ilk)
         kwargs['WD_ilk'] = WD_ilk
@@ -1016,6 +1057,5 @@
                 if any([k in modified_input_dict for k in ['x_ilk', 'y_ilk']]):
-                    z_ilk = self.site.elevation(model_kwargs['x_ilk'], model_kwargs['y_ilk'])
-                    self.site.distance.setup(model_kwargs['x_ilk'], model_kwargs['y_ilk'], model_kwargs['h_ilk'], z_ilk)
-                    model_kwargs.update({k: v for k, v in zip(['dw_ijlk', 'hcw_ijlk', 'dh_ijlk'],
-                                                              self.site.distance(wd_l=wd, WD_ilk=WD_ilk))})
+                    model_kwargs.update({k: v for k, v in zip(
+                        ['dw_ijlk', 'hcw_ijlk', 'dh_ijlk'],
+                        self.site.distance(*[model_kwargs[k + '_ilk'] for k in 'xyh'], wd_l=wd, WD_ilk=WD_ilk))})
                     model_kwargs['cw_ijlk'] = gradients.hypot(model_kwargs['dh_ijlk'], model_kwargs['hcw_ijlk'])
@@ -1046,2 +1086,7 @@ if not self.deflectionModel:
 
+            for i, ewf in enumerate(self.externalWindFarms, I - len(self.externalWindFarms)):
+                deficit_jlk = ewf(i=i, l=np.ones(L, dtype=bool), deficit_jlk=None, dst_xyh_jlk=dst_xyh_jlk,
+                                  **model_kwargs)
+                deficit_iilk = item_assign(deficit_iilk, idx=i, values=deficit_jlk)
+
             # set own deficit to 0
@@ -1048,0 +1093,0 @@ deficit_iilk *= i2i_zero

py_wake/wind_farm_models/minimalistic_wind_farm_model.py

@@ -187,2 +187,3 @@ """
         ws_cutout = ws_cutout or getattr(windTurbines, 'ws_cutout', 25)
+        self.externalWindFarms = []
 
@@ -189,0 +190,0 @@ MinimalisticPredictionModel.__init__(self, correction_factor, latitude, CP=max_cp,

py_wake/wind_farm_models/wind_farm_model.py

@@ -1,17 +0,19 @@
-from abc import abstractmethod, ABC
-from py_wake.site._site import Site, LocalWind
-from py_wake.wind_turbines import WindTurbines
-from py_wake import np
-from py_wake.flow_map import FlowMap, HorizontalGrid, FlowBox, Grid
+import multiprocessing
+from abc import ABC, abstractmethod
+from collections.abc import Iterable
+
 import xarray as xr
-from py_wake.utils import xarray_utils, weibull  # register ilk function @UnusedImport
+from numpy import atleast_1d
 from numpy import newaxis as na
-from numpy import atleast_1d
-from py_wake.utils.model_utils import check_model, fix_shape
-import multiprocessing
-from py_wake.utils.parallelization import get_pool_map, get_pool_starmap, get_map_func
+
+from py_wake import np
+from py_wake.flow_map import FlowBox, FlowMap, Grid, HorizontalGrid
+from py_wake.noise_models.iso import ISONoiseModel
+from py_wake.site._site import LocalWind, Site
+from py_wake.utils import weibull, xarray_utils  # register ilk function @UnusedImport
 from py_wake.utils.functions import arg2ilk
 from py_wake.utils.gradients import autograd
-from py_wake.noise_models.iso import ISONoiseModel
-from collections.abc import Iterable
+from py_wake.utils.model_utils import check_model, fix_shape
+from py_wake.utils.parallelization import get_pool_map, get_pool_starmap
+from py_wake.wind_turbines import WindTurbines
 
@@ -64,2 +66,11 @@
         h, _ = self.windTurbines.get_defaults(len(x), type, h)
+        if len(self.externalWindFarms):
+            err_msg = 'Inflow dependent positions not supported in combination with external wind farms'
+            assert len(np.shape(x)) == 1 or np.shape(x)[1:] == (1, 1), err_msg
+            assert len(np.shape(y)) == 1 or np.shape(y)[1:] == (1, 1), err_msg
+            assert len(np.shape(h)) == 1 or np.shape(h)[1:] == (1, 1), err_msg
+            x = np.r_[x, [wf.wf_x for wf in self.externalWindFarms]]
+            y = np.r_[y, [wf.wf_y for wf in self.externalWindFarms]]
+            h = np.r_[h, [np.mean(wf.windTurbines.hub_height(wf.windTurbines.types()))
+                          for wf in self.externalWindFarms]]
         wd, ws = self.site.get_defaults(wd, ws)
@@ -134,6 +145,5 @@ I, L, K, = len(x), len(np.atleast_1d(wd)), (1, len(np.atleast_1d(ws)))[time is False]
         """
-        if isinstance(time, Iterable) and (
-            len(ws if ws is not None else []) != len(time) or
-            len(wd if wd is not None else []) != len(time)
-        ):  # avoid redundant passing wd & ws for time series sites
+        if isinstance(time, Iterable) and (len(ws if ws is not None else []) != len(time) or
+                                           len(wd if wd is not None else []) != len(time)
+                                           ):  # avoid redundant passing wd & ws for time series sites
             wd = np.zeros(len(time))
@@ -293,3 +303,10 @@ ws = np.zeros(len(time))
         ws_i = np.linspace(0, len(ws) + 1, ws_chunks + 1).astype(int)
+        if n_cpu > 1:
+            map_func = get_pool_map(n_cpu)
+        else:
+            from tqdm import tqdm
 
+            def map_func(f, iter):
+                return tqdm(map(f, iter), total=len(iter), disable=not self.verbose)
+
         if time is False:
@@ -302,4 +319,2 @@ # (wd x ws) matrix
             # (wd, ws) vector
-            if time is True:
-                time = np.arange(len(wd))
             slice_lst = [(slice(wd_i0, wd_i1), slice(wd_i0, wd_i1))
@@ -333,3 +348,3 @@ for wd_i0, wd_i1 in zip(wd_i[:-1], wd_i[1:])]
 
-        return get_map_func(n_cpu, self.verbose), arg_lst, wd_chunks, ws_chunks
+        return map_func, arg_lst, wd_chunks, ws_chunks
 
@@ -668,12 +683,11 @@ def _aep_chunk_wrapper(self, aep_function,
     def flow_box(self, x, y, h, wd=None, ws=None, time=None):
+        if wd is None:
+            wd = self.wd
+        if ws is None:
+            ws = self.ws
         X, Y, H = np.meshgrid(x, y, h)
         x_j, y_j, h_j = X.flatten(), Y.flatten(), H.flatten()
-
-        if time is not None:
-
-            lw_j, WS_eff_jlk, TI_eff_jlk = self.windFarmModel._flow_map(x_j, y_j, h_j, self.sel(time=time))
-        else:
-            wd, ws = self._wd_ws(wd, ws)
-            lw_j, WS_eff_jlk, TI_eff_jlk = self.windFarmModel._flow_map(x_j, y_j, h_j, self.sel(wd=wd, ws=ws))
-
+        wd, ws, sim_res = self._get_flow_map_args(wd, ws, time)
+        lw_j, WS_eff_jlk, TI_eff_jlk = self.windFarmModel._flow_map(x_j[:, na], y_j[:, na], h_j[:, na],
+                                                                    self.localWind, wd, ws, sim_res)
         return FlowBox(self, X, Y, H, lw_j, WS_eff_jlk, TI_eff_jlk)
@@ -704,7 +718,6 @@
                     elif k in {'dw_ijlk', 'hcw_ijlk', 'cw_ijlk', 'dh_ijlk'}:
-                        z_ilk = self.windFarmModel.site.elevation(self.x.ilk(), self.y.ilk())
-                        self.windFarmModel.site.distance.setup(self.x.ilk(), self.y.ilk(), self.h.ilk(), z_ilk)
-                        dist = {k: v for k, v in zip(['dw_ijlk', 'hcw_ijlk', 'cw_ijlk'],
-                                                     self.windFarmModel.site.distance(WD_ilk=self.WD.ilk(),
-                                                                                      wd_l=self.wd.values))}
+                        dist = {k: v for k, v in zip(
+                            ['dw_ijlk', 'hcw_ijlk', 'cw_ijlk'],
+                            self.windFarmModel.site.distance(self.x.ilk(), self.y.ilk(), self.h.ilk(),
+                                                             WD_ilk=self.WD.ilk(), wd_l=self.wd.values))}
                         wt_kwargs.update({k: v for k, v in dist.items() if k in lst})
@@ -724,3 +737,3 @@ if k == 'cw_ijlk':  # pragma: no cover
             setattr(sim_res, k, getattr(self, k))
-        aep_j = self.windFarmModel._aep_map(x_j, y_j, h_j, type, sim_res, n_cpu, memory_GB)
+        aep_j = self.windFarmModel._aep_map(x_j[:, na], y_j[:, na], h_j[:, na], type, sim_res, n_cpu, memory_GB)
         if normalize_probabilities:
@@ -774,4 +787,12 @@ lw_j = self.windFarmModel.site.local_wind(x=x_j, y=y_j, h=h_j, wd=wd, ws=ws)
         X, Y, x_j, y_j, h_j, plane = self._get_grid(grid)
+        wd, ws, sim_res = self._get_flow_map_args(wd, ws, time)
+        lw_j, WS_eff_jlk, TI_eff_jlk = self.windFarmModel._flow_map(
+            x_j[:, na], y_j[:, na], h_j[:, na], self.localWind, wd, ws, sim_res, D_dst=D_dst,
+            memory_GB=memory_GB, n_cpu=n_cpu)
+        return FlowMap(sim_res, X, Y, lw_j, WS_eff_jlk, TI_eff_jlk, plane=plane)
+
+    def _get_flow_map_args(self, wd, ws, time):
         if 'time' in self:
             sim_res = self.sel(time=(time, slice(time))[time is None])
+            wd, ws = sim_res.wd.values, sim_res.ws.values
         else:
@@ -786,7 +807,4 @@ if wd is None:
             setattr(sim_res, k, getattr(self, k))
+        return wd, ws, sim_res
 
-        lw_j, WS_eff_jlk, TI_eff_jlk = self.windFarmModel._flow_map(
-            x_j, y_j, h_j, sim_res, D_dst=D_dst, memory_GB=memory_GB, n_cpu=n_cpu)
-        return FlowMap(sim_res, X, Y, lw_j, WS_eff_jlk, TI_eff_jlk, plane=plane)
-
     def _wd_ws(self, wd, ws):
@@ -836,7 +854,9 @@ if wd is None:
     if __name__ == '__main__':
-        from py_wake.examples.data.iea37 import IEA37Site, IEA37_WindTurbines
-        from py_wake import IEA37SimpleBastankhahGaussian
         import warnings
+
         import matplotlib.pyplot as plt
 
+        from py_wake import IEA37SimpleBastankhahGaussian
+        from py_wake.examples.data.iea37 import IEA37_WindTurbines, IEA37Site
+
         site = IEA37Site(16)
@@ -843,0 +863,0 @@ x, y = site.initial_position.T

py_wake/wind_turbines/_wind_turbines.py

@@ -153,3 +153,3 @@ from py_wake import np
 
-    def plot_xy(self, x, y, types=None, wd=None, yaw=0, tilt=0, normalize_with=1, ax=None):
+    def plot_xy(self, x, y, types=None, wd=None, yaw=0, tilt=0, normalize_with=1, ax=None, wt_number=1):
         """Plot wind farm layout including type name and diameter
@@ -180,3 +180,2 @@
             ax = plt.gca()
-        markers = np.array(list("213v^<>o48spP*hH+xXDd|_"))
 
@@ -188,3 +187,2 @@ assert len(x) == len(y)
         tilt = np.zeros_like(x) + tilt
-
         x, y, D = [np.asarray(v) / normalize_with for v in [x, y, self.diameter(types)]]
@@ -206,3 +204,3 @@ R = D / 2
 
-        for t, m in zip(np.unique(types), markers):
+        for t in np.unique(types):
             color = cmap[t % cmap.shape[0], :]
@@ -212,4 +210,5 @@ # ax.plot(np.asarray(x)[types == t], np.asarray(y)[types == t], '%sk' % m, label=self._names[int(t)])
 
-        for i, (x_, y_, r) in enumerate(zip(x, y, R)):
-            ax.annotate(i, (x_ + r, y_ + r), fontsize=7)
+        if wt_number:
+            for i, (x_, y_, r) in enumerate(zip(x, y, R)):
+                ax.annotate(i, (x_ + r, y_ + r), fontsize=7)
         # ax.legend(loc=1)
@@ -246,4 +245,2 @@
             ax = plt.gca()
-        markers = np.array(list("213v^<>o48spP*hH+xXDd|_"))
-        colors = ['k', 'gray', 'r', 'g', 'k'] * 5
 
@@ -265,3 +262,3 @@ yaw = np.zeros_like(y) + yaw
                 circle = Ellipse((y_, h_ + z_), d * np.sin(np.deg2rad(wd - yaw_)),
-                                 d, angle=-tilt_, ec=colors[t], fc="None", zorder=32)
+                                 d, angle=-tilt_, ec=cmap[t % cmap.shape[0]], fc="None", zorder=32)
                 ax.add_artist(circle)
@@ -272,3 +269,4 @@ else:
 
-        for t, m, c in zip(np.unique(types), markers, colors):
+        for t in np.unique(types):
+            c = cmap[t % cmap.shape[0]]
             ax.plot([], [], '2', color=c, label=self._names[int(t)])
@@ -280,4 +278,4 @@
 
-    def plot(self, x, y, type=None, wd=None, yaw=0, tilt=0, normalize_with=1, ax=None):
-        return self.plot_xy(x, y, type, wd, yaw, tilt, normalize_with, ax)
+    def plot(self, x, y, type=None, wd=None, yaw=0, tilt=0, normalize_with=1, ax=None, wt_number=1):
+        return self.plot_xy(x, y, type, wd, yaw, tilt, normalize_with, ax, wt_number)
 
@@ -284,0 +282,0 @@ def plot_power_ct(self, ax=None, ws=np.linspace(0, 25, 1000), **wt_kwargs):

RECORD

@@ -7,7 +7,7 @@ docs/__init__.py,sha256=47DEQpj8HBSa-_TImW-5JCeuQeRkm5NMpJWZG3hSuFU,0
 py_wake/__init__.py,sha256=-0bNwPpZWQXG-_zE_7hFRoJeG0SRBuDeGM9n5aFWgS0,1188
-py_wake/flow_map.py,sha256=I5InulXHI0dIP-JXuJlnj9jyHWNeOsnpKfVtuVv44IQ,22119
+py_wake/flow_map.py,sha256=zUFPSUQZmxW5LwFTOiBdaL9fPRJxt87UGGemc_PGl3o,22411
 py_wake/superposition_models.py,sha256=mlZSy-T6YiAhdSgAu-reP1YECRsPpgHh8GDCfPAsZ6w,10124
-py_wake/version.py,sha256=Vi-fvuY4jNLmkjfsZxZ_bieoAIk2n0J8Fd96j6hAUrI,706
+py_wake/version.py,sha256=7zrBk63M80mRuOeY6dfWyf-ZQjC6DN9bM1NV1oUDz0Q,706
 py_wake/deficit_models/__init__.py,sha256=w1S4F3FmAgtwqGO90fASYoqxOb6SsgjbasiD3V8_1uk,709
-py_wake/deficit_models/deficit_model.py,sha256=BQ8vJ-Nlw4CDWHdqIkNtbAc52OOwfEi4n6P0Vmg8bLg,10255
+py_wake/deficit_models/deficit_model.py,sha256=iLb9wUO8wDeYAEBAoAlS2W4Q4B16ZSO0L4DhVq3NxYQ,10280
 py_wake/deficit_models/fuga.py,sha256=Jhi52MpvF7-KTk62dZfdLLUf4piCslwIBJTr7N17Hng,13863
@@ -278,13 +278,14 @@ py_wake/deficit_models/gaussian.py,sha256=_ntryBrucuo9EMSLnyJH3dt9mM9SnB0DxP9wzH_w7Zo,36119
 py_wake/noise_models/__init__.py,sha256=47DEQpj8HBSa-_TImW-5JCeuQeRkm5NMpJWZG3hSuFU,0
-py_wake/noise_models/iso.py,sha256=-E6X7pEfOpg6QqTGbLEjLY22NAiJg4Ec_VvWvsVUpR8,12438
+py_wake/noise_models/iso.py,sha256=GYaTGuCnjwSlt8et8zxDM21XQhXd0_FHckEkEHnqJMk,12482
 py_wake/rotor_avg_models/__init__.py,sha256=HK_AO6731z8G3EosHb4bcd0nW7YobeCBtrxhBQQ87aY,357
 py_wake/rotor_avg_models/area_overlap_model.py,sha256=NPEpdtZv9QnoXt4DC_1JSdCETFjNzhinFrJZ1D2lk3k,5237
 py_wake/rotor_avg_models/gaussian_overlap_.02_.02_128_512.nc,sha256=qcFEq0C0z60N2UErbYm-8BLqPALFPwqs0CI552U6RoI,4031394
-py_wake/rotor_avg_models/gaussian_overlap_model.py,sha256=R4bstoVSi6baXaCeivlDs0SCAeQwxrz_yK3FrgJIB0c,5622
+py_wake/rotor_avg_models/gaussian_overlap_model.py,sha256=t8UVAKe7aSPNfyvqEZQxb5QwfxyChu8W6zV147GrWnc,5727
 py_wake/rotor_avg_models/rotor_avg_model.py,sha256=mTkz5shWbc48edQHSoNRjxhKyfWR-Aq7a37-g-2TX7k,10959
 py_wake/site/__init__.py,sha256=YeIRWHuiIo0cRD3NIcaqgoWL6u-NfQkDo3uQsgAxjj4,161
-py_wake/site/_site.py,sha256=bwVmGSVrGkxnv8POZJ225enQDcpOoj1SwcnzlAnrcMs,16367
-py_wake/site/distance.py,sha256=-H-xZ5emerf_eO-1dTQ-crG9xToeQUb7DFzUnaZIfek,11365
-py_wake/site/jit_streamline_distance.py,sha256=a2FVowV8cjUb-U18riZ4k2NskGPP0KSyqL4w621A2X8,4925
+py_wake/site/_site.py,sha256=zvIXS6oPbUffMTXAIq1TXXKhYqePqEFNGAejPtPRxyQ,16256
+py_wake/site/distance.py,sha256=61kvn9ye6ld_vn7rTdJ0GgLjA41uB4e_esnoXDWumzw,8947
+py_wake/site/jit_streamline_distance.py,sha256=GJo2mn_bTjuGX3-gSGkv1dkRfo7hGDVKDmefRWYk5o8,159
 py_wake/site/shear.py,sha256=OGU2eQCuZ7H3hAjAbATK9ZXqIWBFjDazp9dmKumDnSw,4213
+py_wake/site/streamline_distance.py,sha256=Qbjv2GCrPzxFlAp6ev5UhQGPfq2CVyFDzcK4iAHWkP0,4434
 py_wake/site/wasp_grid_site.py,sha256=81zqtUIwCAUHaZdgnNvzazBjEIPYk22Gr6r7IQzmsZk,10190
@@ -295,5 +296,5 @@ py_wake/site/xrsite.py,sha256=-Kr2LMx0dpKI-ke4z3hbRyv-GWulLntjkSU-nY1LWIA,19887
 py_wake/tests/notebook.py,sha256=21RVNXcqRCy6b3GoRcKxR5E0IGV2d9U6BaqqOKZ8f8Y,6804
-py_wake/tests/test_flow_map.py,sha256=8lQzTxsQp4xUKkdRGkmdS80wYbN5QcwsYpOSJ0z_vhM,20290
-py_wake/tests/test_main.py,sha256=aouzNZBMsqrRitG9HfEOCPx8Fe6IH0mM_eIyeQ_TqnM,1853
-py_wake/tests/test_notebooks.py,sha256=805s9rBBr1dWEtapvK2BwrfWUPSWWiA48Dz-cwciObI,1956
+py_wake/tests/test_flow_map.py,sha256=V_dZu9y1K4QaAM7whhKyvTG2O57LmLk5o6y08mbSauw,22448
+py_wake/tests/test_main.py,sha256=pc01ZUXMnFDHZwNh7ulqP5FWZDt9I7-gp0D_TN2-9Bs,1914
+py_wake/tests/test_notebooks.py,sha256=j2mlOAzwMwe6PfEFS2sEl4aQJDxfDQHR8DM6gZq9B-Q,1958
 py_wake/tests/test_superposition_models.py,sha256=dxcgeB_LXn1USv_4QNRuyjL83Cy0qhV1yPdJN6enRj4,8759
@@ -304,7 +305,7 @@ py_wake/tests/test_deficit_models/__init__.py,sha256=47DEQpj8HBSa-_TImW-5JCeuQeRkm5NMpJWZG3hSuFU,0
 py_wake/tests/test_deficit_models/test_blockage_models.py,sha256=guS3PX7YdNfI1Po4JzEhGu6BT6ODSxEBR6fBIMUorko,11856
-py_wake/tests/test_deficit_models/test_deficit_models.py,sha256=vwwj3lSs_4kWCl55-qVs0-3bwdmUT7ZXICpVAuZSfRU,33990
+py_wake/tests/test_deficit_models/test_deficit_models.py,sha256=vIe5H_cX0rZ-V1bD3I4FWgnaeR1qy6NGc-uCefA9qhI,34083
 py_wake/tests/test_deficit_models/test_fuga.py,sha256=zuWx-oAThjYGdubhX9c-dR8-wsnBJMTi7uA6R17_kos,24769
 py_wake/tests/test_deficit_models/test_gcl.py,sha256=i2mIAoNUnGswJ8xXVjjHVw9Imgzhs6thfI8zkgLCt04,1189
 py_wake/tests/test_deficit_models/test_noj.py,sha256=qJVv3q3Lp_QILFn-6xdslWevvR4Dmnx0AWj39XCIdXw,5495
-py_wake/tests/test_deficit_models/test_rans_lut.py,sha256=rfvT0178DfxA52wlscSDL56tJapBVI_BJu1CnWGkCQA,17738
+py_wake/tests/test_deficit_models/test_rans_lut.py,sha256=jA8h7Fp4euRvjQNrfnT-gqLGofZNwjQElKBuJnbcIpE,17130
 py_wake/tests/test_deficit_models/test_selfsimilarity.py,sha256=mI5DUstBbaC36JQ-eCOMrE88jubvZjb3tp_DUQQzU2Y,4446
@@ -338,7 +339,7 @@ py_wake/tests/test_deflection_models/__init__.py,sha256=47DEQpj8HBSa-_TImW-5JCeuQeRkm5NMpJWZG3hSuFU,0
 py_wake/tests/test_sites/__init__.py,sha256=47DEQpj8HBSa-_TImW-5JCeuQeRkm5NMpJWZG3hSuFU,0
-py_wake/tests/test_sites/test_distances.py,sha256=AEB6EwwjE3yo5qkwKcbixWa9xtIByalkFsZPl-ITcAE,16023
+py_wake/tests/test_sites/test_distances.py,sha256=VKox97PB3HmHNH0MO9IbVLiniolY2ffwqz9pNrGiu6M,16338
 py_wake/tests/test_sites/test_iea37_readers.py,sha256=nEoXISd3M7USxJe4PFBI2wUlLTjiL6NEnXMvFWHzwxo,1427
 py_wake/tests/test_sites/test_shear_models.py,sha256=g4zg1QCvl2WQY6f10BiRlBAnyJyaaQRVpeiLnkNT4SU,4053
-py_wake/tests/test_sites/test_site.py,sha256=NZ_gm8LvHwY0PYVeJbXjtUONjSGaoNCZQ8Q8Z938rrU,11267
-py_wake/tests/test_sites/test_wasp_grid_site.py,sha256=8Swf470l8PdI9FRs-yF27S4pP7X08kXRFTZCK1YYNlU,15686
+py_wake/tests/test_sites/test_site.py,sha256=NpE0hX1l17fZmzmrBR5x8w12Hc8Bu5w3NzGG9W60URk,11205
+py_wake/tests/test_sites/test_wasp_grid_site.py,sha256=miNHoOVFt31YqM0Auo97GkAfQ86eyUvc_xkoVO86xEA,15662
 py_wake/tests/test_sites/test_xrsite.py,sha256=3lASkg3Cdv-MQwrXZRuGqQDGONy6pzRdLOurQTneKGI,31239
@@ -353,6 +354,6 @@ py_wake/tests/test_turbulence_models/__init__.py,sha256=47DEQpj8HBSa-_TImW-5JCeuQeRkm5NMpJWZG3hSuFU,0
 py_wake/tests/test_utils/test_gpu_utils.py,sha256=vTVMTLrb9KKr7bHr-AbZZwytV9fhqABq0cn9K4-pt1I,1292
-py_wake/tests/test_utils/test_gradient_utils.py,sha256=phZtdTCg64KgjTvyK4n_HHssPJCYAlJUUHUc9GRWBlU,15810
+py_wake/tests/test_utils/test_gradient_utils.py,sha256=vyacCD9xviDv1YUBmgRoQbG4E1ZiMpSg27anfunw3Kw,16514
 py_wake/tests/test_utils/test_grid_interpolator.py,sha256=XWXlIhCpxoRxlSBM6RUeEmEVfYOe3c2wlhjo9hibyWA,7876
 py_wake/tests/test_utils/test_layouts.py,sha256=9soiopXExIdV1jxSnoG6te_54qbxfRGEd3AdpoGmpNM,1199
-py_wake/tests/test_utils/test_model_gradients.py,sha256=Rni1JhZWB_gOmK-aBNFR7bi8q5VlqKX3ftgA6Nbgmc0,12587
+py_wake/tests/test_utils/test_model_gradients.py,sha256=JHZ-uZfNDW-O8sNTssBmE-EJR_XDSAan1OmCtaSROhE,12638
 py_wake/tests/test_utils/test_model_utils.py,sha256=219lb_0_74xshlw2lKdEgTyPGSRXe8F-uplifqpRdnM,3847
@@ -363,3 +364,3 @@ py_wake/tests/test_utils/test_most.py,sha256=WG2Ng8T1mcmIAsP23OnQeGCoVQ9D8NEX-2fmhH4tMd8,1311
 py_wake/tests/test_utils/test_plotting.py,sha256=tCmEUdvZinzcFMiBGRth8MbFKuVSxQFDuPNh762GBps,286
-py_wake/tests/test_utils/test_profiling.py,sha256=g65cCmFOxs2xsQ9PVJTa6TopSpbRGRI-F3e3zpA10Ts,1562
+py_wake/tests/test_utils/test_profiling.py,sha256=jwanMam7pqEVzlKdsjKqDzWJN1Y7PqY0RN6FqYggNIk,1558
 py_wake/tests/test_utils/test_xarray_utils.py,sha256=abapAtIBQ931ogNqHvi6X2xL_gAKv5DXMdLV1pqiMss,3069
@@ -374,6 +375,7 @@ py_wake/tests/test_utils/test_floris_api/test_floris_turbine.py,sha256=yrLQsNEyV3tiTFVLIxOaNEiAGhYPf3_ClI7Vbsuz9Vc,317
 py_wake/tests/test_wind_farm_models/test_enginering_wind_farm_model.py,sha256=IvbRInNWqnC8FVk8A-IwmDiYOZ3oGNR4Z8drlRHHpnk,28773
+py_wake/tests/test_wind_farm_models/test_external_windfarms.py,sha256=MszAdqFAhnQyxbIXABMc-qgGmZ6y_vsL-XMDleiiKVU,12489
 py_wake/tests/test_wind_farm_models/test_memory_usage.py,sha256=nDGlWHDNVAUrZZD9BOOC-AluRWKuHCLiljp3m8ZofJs,2016
 py_wake/tests/test_wind_farm_models/test_minimalistic_model.py,sha256=Nmu7eNldVD_m_G7bmsOf8eUGkYf_bF5kUhNwFWF9fj8,1566
 py_wake/tests/test_wind_farm_models/test_parallel_aep.py,sha256=QvBk7bFpaM5IGRBhhZmUN64SI1eRlm7XOUx3JecpHwY,3037
-py_wake/tests/test_wind_farm_models/test_parallel_daep.py,sha256=B128GWrFQl2D8a30-Fkn1emqa4gTud2icNo1JVyYtrY,3272
+py_wake/tests/test_wind_farm_models/test_parallel_daep.py,sha256=NrbM7uF6Fed89RqTiItsErEDBB_UuWezYzzLqkK1GKg,3309
 py_wake/tests/test_wind_farm_models/test_simulation_result.py,sha256=VoffAnO6EwA4J_4-AbgoCgrU5Lbj-3f_HbC-YeBVav8,1496
@@ -400,15 +402,15 @@ py_wake/tests/test_wind_farm_models/test_wind_farm_model.py,sha256=hsv11zdQD4X5p1aB09vOiJKJx1ZxlqqzZ1QnqBDoUB8,14464
 py_wake/utils/gpu_utils.py,sha256=qun-tz8HuwLForxFk7H8DlkdPNmjeRiqMew7f6LOc1c,1618
-py_wake/utils/gradients.py,sha256=gb_GoWJWx38vQN0z6HKWB62AQPnDo53oHUDEhcLd16s,14498
+py_wake/utils/gradients.py,sha256=uuhBCbpb7y9K02qbW0m6Ryr-QI7HCbsQzD3Q2uFwbLw,15207
 py_wake/utils/grid_interpolator.py,sha256=Yq6VFkEMkxJkALC74ATvvHYpTlRMsinkG6VTe_c49Pg,8231
 py_wake/utils/ieawind37_utils.py,sha256=t32eUKrPP6oH2G7Lc53l03TZOe7cj4XhBt_FlNZB5kQ,305
 py_wake/utils/layouts.py,sha256=YHV1FwZJECBjTqqYIi8ipVIOefXYoWS7AarK7ND048M,998
-py_wake/utils/model_utils.py,sha256=JYC1xYuJ5U720Tc7_jB_iDFvYaRZ_WlN-Py0JPOho_M,15389
+py_wake/utils/model_utils.py,sha256=CkdZ38hwCkmyfTO5x6iwAgrVi3o0hOohDVZ213ti_Uw,15350
 py_wake/utils/most.py,sha256=9DEZ5C06qdttfGdfUCyaFCGt-YT8xPnQ3URIuHpbcdE,1208
 py_wake/utils/notebook_extensions.py,sha256=E0UOiQ1AAMaIbnnzsyj-gGv8oRWD8wAvHsj2BN10LXU,704
 py_wake/utils/numpy_utils.py,sha256=I9URsvF0ALXOFwdlkTK-mQLul0d-5wLjSWuEeAmE5ZU,3192
-py_wake/utils/parallelization.py,sha256=7eCsS-2A26N8P8dmdkX2d7eW4jy0jQoFOcPxA24DbcU,2082
+py_wake/utils/parallelization.py,sha256=k_w4AQMndE4acS0vwBMlZx4AgT_Oi7N7OK3OqVJzc5I,2101
 py_wake/utils/plotting.py,sha256=zXharIIWIwl-E0Xh2JP3qAnYEEAKGJidQcldhRd0dRo,568
 py_wake/utils/profiling.py,sha256=GuE0UcQaZ3KsVnklLHSDFpCc0jQ4fixh0VREwA3aHQ0,4618
-py_wake/utils/rans_lut_utils.py,sha256=nQut0o1agTcP-vAumQ9DjITtO5oJUD0NB4S0nVhApoQ,11381
-py_wake/utils/streamline.py,sha256=NSiA7acgkeh8iUepqalvSt7XhWS7U4POvoZq2LZCaIc,2172
+py_wake/utils/rans_lut_utils.py,sha256=-zs25tCdg2heaoCjyEDdk0vSDTImdhjtTpbDH_YEClI,12423
+py_wake/utils/streamline.py,sha256=hFIvkHsdUPLzySvQYdpX6DEUMaz5NayQBEg2fuppPh4,2574
 py_wake/utils/tensorflow_surrogate_utils.py,sha256=QdIPqxoYiuXFvieHQyP-TE6gS8syXbxoUZ1liCBFil0,7198
@@ -423,7 +425,8 @@ py_wake/utils/weibull.py,sha256=cD1SIyXHEn17IQueZnsy9-upo2oZivVCpYOWka-_DsE,1262
 py_wake/wind_farm_models/__init__.py,sha256=cl0nC9vqhXgpRuaCFkQpDET2Yyzbwnd5Kv2wgZHUI6E,137
-py_wake/wind_farm_models/engineering_models.py,sha256=V9Ebkm2uMz0mo4iyzu_ELT4pEMn7qs-uRsJi-q_esbI,60738
-py_wake/wind_farm_models/minimalistic_wind_farm_model.py,sha256=geyeFlxI71PPs19pDSWDaQgiSyxaJqNTcuXI1he93_0,11124
-py_wake/wind_farm_models/wind_farm_model.py,sha256=7_Bcwlca0Q-EmoFcsM8vFtlQhV-dEcByTtY5LLUlsA4,41812
+py_wake/wind_farm_models/engineering_models.py,sha256=heyGqJroRPO_BjRaMf1QeEyu-kpiUkxntQTw4r5VfeA,63601
+py_wake/wind_farm_models/external_wind_farm_models.py,sha256=8Qjkzf9w9RFSVr7AHEE3wPtgiwlsZue4SGlBhb1w32k,12649
+py_wake/wind_farm_models/minimalistic_wind_farm_model.py,sha256=9BepPjAo_WDMcG4feUQo7YDY2WSFJOpNMmzgT0yNBn4,11160
+py_wake/wind_farm_models/wind_farm_model.py,sha256=QvYmeqAVz2QyrjH21VEt1T52lJS2NQcTq3cW-W1qUFQ,42814
 py_wake/wind_turbines/__init__.py,sha256=w1D9rLfxk7m_UdrqbVVokWwAikxIzeRh6Wb9zVT2Mhs,145
-py_wake/wind_turbines/_wind_turbines.py,sha256=IAIaERtPtRR6QwhHweKC1Vqi8IYqjBV2ebmaby2Mqqo,19010
+py_wake/wind_turbines/_wind_turbines.py,sha256=ukh8Kw4tUHH5vkXiVC8JaI-5xuIjCRVUROWInnJQJcw,18915
 py_wake/wind_turbines/generic_wind_turbines.py,sha256=SiV-hganUn9TO8ItM8ZPVuJ0BqKIuig5-vA6miSGJeY,9456
@@ -433,6 +436,6 @@ py_wake/wind_turbines/power_ct_functions.py,sha256=CjWKehvop9kMPclU54pFp_hRG7_3bqCx16cJlz5QduU,22663
 py_wake/wind_turbines/wind_turbines_deprecated.py,sha256=HpNmBR8CJL4-8JBaygDI0t086qfw5bR2DOQI8Ox4AZ4,6250
-py_wake-2.6.13.dist-info/licenses/LICENSE,sha256=XE2CGPqQgzSXqIajXpAVYJ5SRNmaWOIeMePK6MocsuY,1084
-py_wake-2.6.13.dist-info/METADATA,sha256=fpTWS7ZgV4KrfCqxIDfNfV-qTxPc295iQmBlADMgfpg,3596
-py_wake-2.6.13.dist-info/WHEEL,sha256=_zCd3N1l69ArxyTb8rzEoP9TpbYXkqRFSNOD5OuxnTs,91
-py_wake-2.6.13.dist-info/top_level.txt,sha256=GsaXU4YwyMkZZ6dkb4h0FMc5RaLIT2Qns_YoScKoXdk,20
-py_wake-2.6.13.dist-info/RECORD,,
+py_wake-2.6.14.dist-info/licenses/LICENSE,sha256=XE2CGPqQgzSXqIajXpAVYJ5SRNmaWOIeMePK6MocsuY,1084
+py_wake-2.6.14.dist-info/METADATA,sha256=GTsCfUqg3V5gNXNjdrnN3JYWPDYHOun_fqN0ZCqDdj0,3596
+py_wake-2.6.14.dist-info/WHEEL,sha256=_zCd3N1l69ArxyTb8rzEoP9TpbYXkqRFSNOD5OuxnTs,91
+py_wake-2.6.14.dist-info/top_level.txt,sha256=GsaXU4YwyMkZZ6dkb4h0FMc5RaLIT2Qns_YoScKoXdk,20
+py_wake-2.6.14.dist-info/RECORD,,

Alert delta unavailable

Currently unable to show alert delta for PyPI packages.

No metric changes