engineering_models.py

from abc import abstractmethod
from numpy import newaxis as na
import numpy as np
from py_wake.deficit_models import DeficitModel
from py_wake.superposition_models import SuperpositionModel
from py_wake.turbulence_models.turbulence_model import TurbulenceModel
from py_wake.wind_farm_models.wind_farm_model import WindFarmModel
from py_wake.deflection_models.deflection_model import DeflectionModel
from py_wake.gradients import use_autograd_in, autograd


class EngineeringWindFarmModel(WindFarmModel):
    """
    Base class for engineering wake models

    General suffixes:

    - i: turbines ordered by id
    - j: downstream points/turbines
    - k: wind speeds
    - l: wind directions

    Arguments available for calc_deficit (specifiy in args4deficit):

    - WS_ilk: Local wind speed without wake effects
    - TI_ilk: local turbulence intensity without wake effects
    - WS_eff_ilk: Local wind speed with wake effects
    - TI_eff_ilk: local turbulence intensity with wake effects
    - D_src_il: Diameter of source turbine
    - D_dst_ijl: Diameter of destination turbine
    - dw_ijlk: Downwind distance from turbine i to point/turbine j
    - hcw_ijlk: Horizontal cross wind distance from turbine i to point/turbine j
    - dh_ijl: vertical distance from turbine i to point/turbine j
    - cw_ijlk: Cross wind(horizontal and vertical) distance from turbine i to point/turbine j
    - ct_ilk: Thrust coefficient

    """
    default_grid_resolution = 500

    def __init__(self, site, windTurbines, wake_deficitModel, superpositionModel,
                 blockage_deficitModel=None, deflectionModel=None, turbulenceModel=None):

        WindFarmModel.__init__(self, site, windTurbines)

        assert isinstance(wake_deficitModel, DeficitModel)
        assert isinstance(superpositionModel, SuperpositionModel)
        assert blockage_deficitModel is None or isinstance(blockage_deficitModel, DeficitModel)
        assert deflectionModel is None or isinstance(deflectionModel, DeflectionModel)
        assert turbulenceModel is None or isinstance(turbulenceModel, TurbulenceModel)
        self.site = site
        self.windTurbines = windTurbines
        self.wake_deficitModel = wake_deficitModel
        self.superpositionModel = superpositionModel
        self.blockage_deficitModel = blockage_deficitModel
        self.deflectionModel = deflectionModel
        self.turbulenceModel = turbulenceModel

        self.wec = 1  # wake expansion continuation (wake-width scale factor) see
        # Thomas, J. J. and Ning, A., “A Method for Reducing Multi-Modality in the Wind Farm Layout Optimization Problem,”
        # Journal of Physics: Conference Series, Vol. 1037, The Science of Making
        # Torque from Wind, Milano, Italy, jun 2018, p. 10.
        self.deficit_initalized = False

        self.args4deficit = self.wake_deficitModel.args4deficit
        if self.blockage_deficitModel:
            self.args4deficit = set(self.args4deficit) | set(self.blockage_deficitModel.args4deficit)

    def __str__(self):
        def name(o):
            return o.__class__.__name__

        models = [self.__class__.__bases__[0].__name__,
                  "%s-wake" % name(self.wake_deficitModel)]
        if self.blockage_deficitModel:
            models.append("%s-blockage" % name(self.blockage_deficitModel))
        models.append("%s-superposition" % (name(self.superpositionModel)))
        if self.deflectionModel:
            models.append("%s-deflection" % name(self.deflectionModel))
        if self.turbulenceModel:
            models.append("%s-turbulence" % name(self.turbulenceModel))
        return "%s(%s)" % (name(self), ", ".join(models))

    def _init_deficit(self, **kwargs):
        """Calculate layout dependent wake (and blockage) deficit terms"""
        self.wake_deficitModel._calc_layout_terms(**kwargs)
        self.wake_deficitModel.deficit_initalized = True
        if self.blockage_deficitModel:
            self.blockage_deficitModel._calc_layout_terms(**kwargs)
            self.blockage_deficitModel.deficit_initalized = True

    def _reset_deficit(self):
        self.wake_deficitModel.deficit_initalized = False
        if self.blockage_deficitModel:
            self.blockage_deficitModel.deficit_initalized = False

    def _calc_deficit(self, dw_ijlk, **kwargs):
        """Calculate wake (and blockage) deficit"""
        deficit = self.wake_deficitModel.calc_deficit(dw_ijlk=dw_ijlk, **kwargs)

        # the split line between wake and blockage is set slightly downstream to handle
        # numerical inaccuracy in the trigonometric functions that calculates dw_ijlk
        rotor_pos = 1e-10
        if self.blockage_deficitModel is None:
            deficit *= (dw_ijlk > 0)
        elif self.blockage_deficitModel != self:
            # downstream wake deficit + upstream blockage
            deficit = ((dw_ijlk > rotor_pos) * deficit +
                       (dw_ijlk <= rotor_pos) * self.blockage_deficitModel.calc_deficit(dw_ijlk=dw_ijlk, **kwargs))

        return deficit

    def calc_wt_interaction(self, x_i, y_i, h_i=None, type_i=0, wd=None, ws=None, yaw_ilk=None):
        """See WindFarmModel.calc_wt_interaction"""
        type_i, h_i, D_i = self.windTurbines.get_defaults(len(x_i), type_i, h_i)
        wd, ws = self.site.get_defaults(wd, ws)

        # Find local wind speed, wind direction, turbulence intensity and probability
        lw = self.site.local_wind(x_i=x_i, y_i=y_i, h_i=h_i, wd=wd, ws=ws)

        # Calculate down-wind and cross-wind distances
        dw_iil, hcw_iil, dh_iil, dw_order_indices_dl = self.site.wt2wt_distances(x_i, y_i, h_i, lw.WD_ilk.mean(2))
        self._validate_input(dw_iil, hcw_iil)

        I, L = dw_iil.shape[1:]
        K = lw.WS_ilk.shape[2]
        WS_eff_ilk = lw.WS_ilk.copy()
        TI_eff_ilk = lw.TI_ilk.copy()
        if yaw_ilk is None:
            yaw_ilk = np.zeros((I, L, K))
        else:
            yaw_ilk = np.deg2rad(yaw_ilk)

        if self.wec != 1:
            hcw_iil = hcw_iil / self.wec

        # add eps to avoid non-differentiable 0

        cw_iil = np.sqrt(hcw_iil**2 + dh_iil**2)

        kwargs = {'localWind': lw,
                  'WS_eff_ilk': WS_eff_ilk, 'TI_eff_ilk': TI_eff_ilk,
                  'type_i': type_i, 'h_i': h_i, 'D_i': D_i, 'yaw_ilk': yaw_ilk,
                  'dw_iil': dw_iil, 'hcw_iil': hcw_iil, 'cw_iil': cw_iil, 'dh_iil': dh_iil,
                  'dw_order_indices_dl': dw_order_indices_dl, 'I': I, 'L': L, 'K': K}
        WS_eff_ilk, TI_eff_ilk, power_ilk, ct_ilk = self._calc_wt_interaction(**kwargs)
        return WS_eff_ilk, TI_eff_ilk, power_ilk, ct_ilk, lw

    @abstractmethod
    def _calc_wt_interaction(self, **kwargs):
        """calculate WT interaction"""

    def _flow_map(self, x_j, y_j, h_j,
                  wt_x_i, wt_y_i, wt_h_i, wt_type_i, yaw_ilk,
                  WD_ilk, WS_ilk, TI_ilk, WS_eff_ilk, TI_eff_ilk, ct_ilk,
                  wd, ws):
        """call this function via SimulationResult.flow_map"""
        # calculate distances
        wt_d_i = self.windTurbines.diameter(wt_type_i)

        lw_j = self.site.local_wind(x_i=x_j, y_i=y_j, h_i=h_j, wd=wd, ws=ws)
        if len(wt_x_i) == 0:
            # If not turbines just return local wind
            return lw_j, lw_j.WS_ilk, lw_j.TI_ilk

        I, J, L, K = [len(x) for x in [wt_x_i, x_j, wd, ws]]
        WS_eff_jlk = np.zeros((len(x_j), L, K))
        TI_eff_jlk = np.zeros((len(x_j), L, K))

        if self.deflectionModel:
            if yaw_ilk is None:
                yaw_ilk = np.zeros((I, L, K))
            yaw_ilk = np.deg2rad(yaw_ilk)
        if L > 1:
            print("|" + "-" * 100 + "|\n|", end="", flush=True)
        for l in range(L):
            if L > 1 and (l * 100) // L > ((l - 1) * 100) // L:
                print(".", end="", flush=True)
            dw_ijl, hcw_ijl, dh_ijl, _ = self.site.distances(wt_x_i, wt_y_i, wt_h_i, x_j, y_j, h_j,
                                                             wd_il=WD_ilk[:, l:l + 1, :].mean(2))

            if self.wec != 1:
                hcw_ijl = hcw_ijl / self.wec

            if I * J * K * 8 / 1024**2 > 10:
                # one wt at the time to avoid memory problems
                deficit_ijk = np.zeros((I, J, K))
                add_turb_ijk = np.zeros((I, J, K))
                for i in range(I):
                    arg_funcs = {'WS_ilk': lambda: WS_ilk[i, l][na, na],
                                 'WS_eff_ilk': lambda: WS_eff_ilk[i, l][na, na],
                                 'TI_ilk': lambda: TI_ilk[i, l][na, na],
                                 'TI_eff_ilk': lambda: TI_eff_ilk[i, l][na, na],
                                 'yaw_ilk': lambda: yaw_ilk[i, l][na, na],
                                 'D_src_il': lambda: wt_d_i[i][na, na],
                                 'D_dst_ijl': lambda: None,
                                 'dh_ijl': lambda: dh_ijl[i][na],
                                 'h_il': lambda: wt_h_i[i][na, na],
                                 'ct_ilk': lambda: ct_ilk[i, l][na, na]}

                    if self.deflectionModel:
                        dw_ijlk, hcw_ijlk = self.deflectionModel.calc_deflection(
                            dw_ijl[i][na], hcw_ijl[i][na],
                            **{k: arg_funcs[k]() for k in self.deflectionModel.args4deflection})
                    else:
                        dw_ijlk, hcw_ijlk = dw_ijl[i][na, :, :, na], hcw_ijl[i][na, :, :, na]
                    arg_funcs['cw_ijlk'] = lambda: np.hypot(dh_ijl[i][na, :, :, na], hcw_ijlk)
                    arg_funcs['hcw_ijlk'] = lambda: hcw_ijlk

                    args = {k: arg_funcs[k]() for k in self.args4deficit if k != 'dw_ijlk'}
                    deficit_ijk[i] = self._calc_deficit(dw_ijlk=dw_ijlk, **args)[0, :, 0]
                    if self.turbulenceModel:
                        arg_funcs['wake_radius_ijlk'] = lambda: self.wake_deficitModel.wake_radius(
                            dw_ijlk=dw_ijlk, **args)
                        args.update({k: arg_funcs[k]() for k in self.turbulenceModel.args4addturb
                                     if k not in self.args4deficit})
                        add_turb_ijk[i] = self.turbulenceModel.calc_added_turbulence(dw_ijlk=dw_ijlk, **args)[0, :, 0]

            else:

                arg_funcs = {'WS_ilk': lambda: WS_ilk[:, l][:, na],
                             'WS_eff_ilk': lambda: WS_eff_ilk[:, l][:, na],
                             'TI_ilk': lambda: TI_ilk[:, l][:, na],
                             'TI_eff_ilk': lambda: TI_eff_ilk[:, l][:, na],
                             'yaw_ilk': lambda: yaw_ilk[:, l][:, na],
                             'D_src_il': lambda: wt_d_i[:, na],
                             'D_dst_ijl': lambda: None,
                             'dh_ijl': lambda: dh_ijl,
                             'h_il': lambda: wt_h_i[:, na],
                             'ct_ilk': lambda: ct_ilk[:, l][:, na]}

                if self.deflectionModel:
                    dw_ijlk, hcw_ijlk = self.deflectionModel.calc_deflection(
                        dw_ijl, hcw_ijl,
                        **{k: arg_funcs[k]() for k in self.deflectionModel.args4deflection})
                else:
                    dw_ijlk, hcw_ijlk = dw_ijl[..., na], hcw_ijl[..., na]
                arg_funcs['cw_ijlk'] = lambda: np.hypot(dh_ijl[..., na], hcw_ijlk)
                arg_funcs['dw_ijlk'] = lambda: dw_ijlk
                arg_funcs['hcw_ijlk'] = lambda: hcw_ijlk

                args = {k: arg_funcs[k]() for k in self.args4deficit if k != 'dw_ijlk'}
                deficit_ijk = self._calc_deficit(dw_ijlk=dw_ijlk, **args)[:, :, 0]
                if self.turbulenceModel:
                    arg_funcs['wake_radius_ijlk'] = lambda: self.wake_deficitModel.wake_radius(
                        dw_ijlk=dw_ijlk, **args)
                    args.update({k: arg_funcs[k]() for k in self.turbulenceModel.args4addturb
                                 if k not in self.args4deficit})
                    add_turb_ijk = self.turbulenceModel.calc_added_turbulence(dw_ijlk=dw_ijlk, **args)[:, :, 0]
            WS_eff_jlk[:, l] = self.superpositionModel.calc_effective_WS(lw_j.WS_ilk[:, l], deficit_ijk)
            if self.turbulenceModel:
                TI_eff_jlk[:, l] = self.turbulenceModel.calc_effective_TI(lw_j.TI_ilk[:, l], add_turb_ijk)
        return lw_j, WS_eff_jlk, TI_eff_jlk

    def _validate_input(self, dw_iil, hcw_iil):
        I_ = dw_iil.shape[0]
        i1, i2, _ = np.where((np.abs(dw_iil) + np.abs(hcw_iil) + np.eye(I_)[:, :, na]) == 0)
        if len(i1):
            msg = "\n".join(["Turbines %d and %d are at the same position" %
                             (i1[i], i2[i]) for i in range(len(i1))])
            raise ValueError(msg)

    def dAEPdn(self, argnum, gradient_method):
        def aep(x, y, h=None, type=0, wd=None, ws=None, yaw_ilk=None):  # @ReservedAssignment
            if gradient_method == autograd:
                with use_autograd_in():
                    return self.__call__(x, y, h, type, wd, ws, yaw_ilk).aep()
            else:
                return self.__call__(x, y, h, type, wd, ws, yaw_ilk).aep()
        return gradient_method(aep, argnum)

    def dAEPdxy(self, gradient_method, normalize_probabilities=False, with_wake_loss=True, gradient_method_kwargs={}):

        def wrap(x, y, h=None, type=0, wd=None, ws=None, yaw_ilk=None):  # @ReservedAssignment
            def aep(x, y, h, type, wd, ws, yaw_ilk):  # @ReservedAssignment
                if gradient_method == autograd:
                    with use_autograd_in():
                        sr = self.__call__(x, y, h, type, wd, ws, yaw_ilk)
                        return sr.aep(normalize_probabilities=normalize_probabilities, with_wake_loss=with_wake_loss)
                else:
                    return self.__call__(x, y, h, type, wd, ws, yaw_ilk).aep(
                        normalize_probabilities=normalize_probabilities, with_wake_loss=with_wake_loss)
            return (gradient_method(aep, 0, **gradient_method_kwargs)(x, y, h, type, wd, ws, yaw_ilk),
                    gradient_method(aep, 1, **gradient_method_kwargs)(x, y, h, type, wd, ws, yaw_ilk))
        return wrap


class PropagateDownwind(EngineeringWindFarmModel):
    """Downstream wake deficits calculated and propagated in downstream direction.
    Very fast, but ignoring blockage effects
    """

    def __init__(self, site, windTurbines, wake_deficitModel, superpositionModel,
                 deflectionModel=None, turbulenceModel=None):
        """Initialize flow model

        Parameters
        ----------
        site : Site
            Site object
        windTurbines : WindTurbines
            WindTurbines object representing the wake generating wind turbines
        wake_deficitModel : DeficitModel
            Model describing the wake(downstream) deficit
        superpositionModel : SuperpositionModel
            Model defining how deficits sum up
        deflectionModel : DeflectionModel
            Model describing the deflection of the wake due to yaw misalignment, sheared inflow, etc.
        turbulenceModel : TurbulenceModel
            Model describing the amount of added turbulence in the wake
        """
        EngineeringWindFarmModel.__init__(self, site, windTurbines, wake_deficitModel, superpositionModel,
                                          blockage_deficitModel=None, deflectionModel=deflectionModel, turbulenceModel=turbulenceModel)

    def _calc_wt_interaction(self, localWind,
                             WS_eff_ilk, TI_eff_ilk,
                             type_i, h_i, D_i, yaw_ilk,
                             dw_iil, hcw_iil, cw_iil, dh_iil, dw_order_indices_dl, I, L, K):
        """
        Additional suffixes:

        - m: turbines and wind directions (il.flatten())
        - n: from_turbines, to_turbines and wind directions (iil.flatten())

        """
        lw = localWind
        deficit_nk = []

        def ilk2mk(x_ilk):
            return x_ilk.astype(np.float).reshape((I * L, K))

        indices = np.arange(I * I * L).reshape((I, I, L))
        TI_mk = ilk2mk(lw.TI_ilk)
        WS_mk = ilk2mk(lw.WS_ilk)
        WS_eff_mk = []
        yaw_mk = ilk2mk(yaw_ilk)
        power_jlk = []
        ct_jlk = []

        if not self.deflectionModel:
            dw_ijlk, hcw_ijlk = dw_iil[..., na], hcw_iil[..., na]

        if self.turbulenceModel:
            add_turb_nk = np.zeros((I * I * L, K))
            TI_eff_mk = ilk2mk(TI_eff_ilk)

        dw_n, hcw_n, cw_n, dh_n = [a.flatten() for a in [dw_iil, hcw_iil, cw_iil, dh_iil]]

        i_wd_l = np.arange(L)

        # Iterate over turbines in down wind order
        for j in range(I):
            i_wt_l = dw_order_indices_dl[:, j]
            m = i_wt_l * L + i_wd_l  # current wt (j'th most upstream wts for all wdirs)

            # generate indexes of up wind(n_uw) and down wind(n_dw) turbines
            n_uw = indices[:, i_wt_l, i_wd_l][dw_order_indices_dl[:, :j].T, np.arange(L)]
            n_dw = indices[i_wt_l, :, i_wd_l][np.arange(L), dw_order_indices_dl[:, j + 1:].T]

            # Calculate effectiv wind speed at current turbines(all wind directions and wind speeds) and
            # look up power and thrust coefficient
            if j == 0:  # Most upstream turbines (no wake)
                WS_eff_lk = WS_mk[m]
                WS_eff_mk.append(WS_eff_lk)
            else:  # 2..n most upstream turbines (wake)
                deficit2WT = np.array([d_nk2[i] for d_nk2, i in zip(deficit_nk, range(j)[::-1])])
                WS_eff_lk = self.superpositionModel.calc_effective_WS(WS_mk[m], deficit2WT)
                WS_eff_mk.append(WS_eff_lk)
                if self.turbulenceModel:
                    TI_eff_mk[m] = self.turbulenceModel.calc_effective_TI(TI_mk[m], add_turb_nk[n_uw])

            ct_lk, power_lk = self.windTurbines._ct_power(WS_eff_lk, type_i[i_wt_l])

            power_jlk.append(power_lk)
            ct_jlk.append(ct_lk)

            if j < I - 1:

                # Calculate required args4deficit parameters
                arg_funcs = {'WS_ilk': lambda: WS_mk[m][na],
                             'WS_eff_ilk': lambda: WS_eff_mk[-1][na],
                             'TI_ilk': lambda: TI_mk[m][na],
                             'TI_eff_ilk': lambda: TI_eff_mk[m][na],
                             'D_src_il': lambda: D_i[i_wt_l][na],
                             'yaw_ilk': lambda: yaw_mk[m][na],
                             'D_dst_ijl': lambda: D_i[dw_order_indices_dl[:, j + 1:]].T[na],
                             'dh_ijl': lambda: dh_n[n_dw][na],
                             'h_il': lambda: h_i[i_wt_l][na],
                             'ct_ilk': lambda: ct_lk[na],
                             'wake_radius_ijlk': lambda: wake_radius_ijlk
                             }

                if self.deflectionModel:
                    dw_ijlk, hcw_ijlk = self.deflectionModel.calc_deflection(
                        dw_ijl=dw_n[n_dw][na], hcw_ijl=hcw_n[n_dw][na],
                        **{k: arg_funcs[k]() for k in self.deflectionModel.args4deflection})

                else:
                    dw_ijlk, hcw_ijlk = dw_n[n_dw][na, :, :, na], hcw_n[n_dw][na, :, :, na],

                arg_funcs['hcw_ijlk'] = lambda: hcw_ijlk
                # sqrt(a**2+b**2) as hypot does not support complex numbers
                arg_funcs['cw_ijlk'] = lambda: np.sqrt(dh_n[n_dw][na, :, :, na]**2 + hcw_ijlk**2)
                args = {k: arg_funcs[k]() for k in self.args4deficit if k != "dw_ijlk"}

                # Calculate deficit

                deficit = self.wake_deficitModel.calc_deficit(dw_ijlk=dw_ijlk, **args)[0]
                deficit_nk.append(deficit)

                if self.turbulenceModel:

                    if 'wake_radius_ijlk' in self.turbulenceModel.args4addturb:
                        wake_radius_ijlk = self.wake_deficitModel.wake_radius(dw_ijlk=dw_ijlk, **args)
                        arg_funcs['wake_radius_ijlk'] = lambda: wake_radius_ijlk

                    turb_args = {k: arg_funcs[k]() for k in self.turbulenceModel.args4addturb if k != "dw_ijlk"}
                    # Calculate added turbulence
                    add_turb_nk[n_dw] = self.turbulenceModel.calc_added_turbulence(dw_ijlk=dw_ijlk, **turb_args)

        WS_eff_jlk, power_jlk, ct_jlk = np.array(WS_eff_mk), np.array(power_jlk), np.array(ct_jlk)

        dw_inv_indices = (np.argsort(dw_order_indices_dl, 1).T * L + np.arange(L)[na]).flatten()
        WS_eff_ilk = WS_eff_jlk.reshape((I * L, K))[dw_inv_indices].reshape((I, L, K))

        power_ilk = power_jlk.reshape((I * L, K))[dw_inv_indices].reshape((I, L, K))
        ct_ilk = ct_jlk.reshape((I * L, K))[dw_inv_indices].reshape((I, L, K))
        if self.turbulenceModel:
            TI_eff_ilk = TI_eff_mk.reshape((I, L, K))

        return WS_eff_ilk, TI_eff_ilk, power_ilk, ct_ilk


class All2AllIterative(EngineeringWindFarmModel):
    """Wake and blockage deficits calculated from all wt to all points of interest (wt/map points).
    The calculations are iteratively repeated until convergence (change of effective wind speed < convergence_tolerance)"""

    def __init__(self, site, windTurbines, wake_deficitModel, superpositionModel,
                 blockage_deficitModel=None, deflectionModel=None, turbulenceModel=None, convergence_tolerance=1e-6):
        """Initialize flow model

        Parameters
        ----------
        site : Site
            Site object
        windTurbines : WindTurbines
            WindTurbines object representing the wake generating wind turbines
        wake_deficitModel : DeficitModel
            Model describing the wake(downstream) deficit
        superpositionModel : SuperpositionModel
            Model defining how deficits sum up
        blockage_deficitModel : DeficitModel
            Model describing the blockage(upstream) deficit
        deflectionModel : DeflectionModel
            Model describing the deflection of the wake due to yaw misalignment, sheared inflow, etc.
        turbulenceModel : TurbulenceModel
            Model describing the amount of added turbulence in the wake
        convergence_tolerance : float
            maximum accepted change in WS_eff_ilk [m/s]
        """
        EngineeringWindFarmModel.__init__(self, site, windTurbines, wake_deficitModel, superpositionModel,
                                          blockage_deficitModel=blockage_deficitModel, deflectionModel=deflectionModel, turbulenceModel=turbulenceModel)
        self.convergence_tolerance = convergence_tolerance

    def _calc_wt_interaction(self, localWind,
                             WS_eff_ilk, TI_eff_ilk,
                             type_i, h_i, D_i, yaw_ilk,
                             dw_iil, hcw_iil, cw_iil, dh_iil, dw_order_indices_dl, I, L, K):
        lw = localWind
        power_ilk = np.zeros((I, L, K))
        WS_eff_ilk_last = WS_eff_ilk.copy()

        ct_ilk = self.windTurbines.ct(lw.WS_ilk, type_i)
        D_src_il = D_i[:, na]
        args = {'WS_ilk': lw.WS_ilk,
                'TI_ilk': lw.TI_ilk,
                'TI_eff_ilk': lw.TI_ilk,
                'yaw_ilk': yaw_ilk,
                'D_src_il': D_src_il,
                'D_dst_ijl': D_src_il[na],
                'cw_ijlk': cw_iil[..., na],
                'dh_ijl': dh_iil,
                'h_il': h_i[:, na]
                }

        # Iterate until convergence
        for j in range(I):

            ct_ilk, power_ilk = self.windTurbines._ct_power(WS_eff_ilk, type_i)
            args['ct_ilk'] = ct_ilk
            args['WS_eff_ilk'] = WS_eff_ilk
            if self.deflectionModel:
                dw_ijlk, hcw_ijlk = self.deflectionModel.calc_deflection(dw_ijl=dw_iil, hcw_ijl=dw_iil, **args)
                args['dw_ijlk'] = dw_ijlk
                args['hcw_ijlk'] = hcw_ijlk
                args['cw_ijlk'] = np.hypot(dh_iil[..., na], hcw_ijlk)
            else:
                args['dw_ijlk'] = dw_iil[..., na]
                args['hcw_ijlk'] = hcw_iil[..., na]
                self._init_deficit(**args)
            if self.turbulenceModel:
                args['TI_eff_ilk'] = TI_eff_ilk
                if 'wake_radius_ijlk' in self.turbulenceModel.args4addturb:
                    args['wake_radius_ijlk'] = self.wake_deficitModel.wake_radius(**args)

            # Calculate deficit
            deficit_iilk = self._calc_deficit(**args)

            # Calculate effective wind speed
            WS_eff_ilk = self.superpositionModel.calc_effective_WS(lw.WS_ilk, deficit_iilk)
            if self.turbulenceModel:
                add_turb_ijlk = self.turbulenceModel.calc_added_turbulence(**args)
                TI_eff_ilk = self.turbulenceModel.calc_effective_TI(lw.TI_ilk, add_turb_ijlk)

            # Check if converged
            diff = np.abs(WS_eff_ilk_last - WS_eff_ilk)
            max_diff = np.max(diff)
            if self.convergence_tolerance and max_diff < self.convergence_tolerance:
                # print("All2AllIterative converge after %d iterations" % j)
                break
            # i_, l_, k_ = list(zip(*np.where(diff == max_diff)))[0]
            # print("Iteration: %d, max diff: %f, WT: %d, WD: %d, WS: %d" % (j, max_diff, i_, l_, WS_ilk[i_, l_, k_]))
            WS_eff_ilk_last = WS_eff_ilk.copy()
        self._reset_deficit()
        return WS_eff_ilk, TI_eff_ilk, power_ilk, ct_ilk


def main():
    if __name__ == '__main__':
        from py_wake.examples.data.iea37 import IEA37Site, IEA37_WindTurbines
        from py_wake.deficit_models.selfsimilarity import SelfSimilarityDeficit

        import matplotlib.pyplot as plt

        site = IEA37Site(16)
        x, y = site.initial_position.T

        windTurbines = IEA37_WindTurbines()
        from py_wake.deficit_models.noj import NOJDeficit
        from py_wake.superposition_models import SquaredSum

        # NOJ wake model
        noj = PropagateDownwind(site, windTurbines, wake_deficitModel=NOJDeficit(), superpositionModel=SquaredSum())

        # NOJ wake and selfsimilarity blockage
        noj_ss = All2AllIterative(site, windTurbines, wake_deficitModel=NOJDeficit(), superpositionModel=SquaredSum(),
                                  blockage_deficitModel=SelfSimilarityDeficit())

        for wm in [noj, noj_ss]:
            plt.figure()
            wm(x=x, y=y, wd=[30], ws=[9]).flow_map().plot_wake_map()
        plt.show()


main()