From 3ee7f5bf82033c72e449bf754278f06fecbf20ea Mon Sep 17 00:00:00 2001
From: "Mads M. Pedersen" <mmpe@dtu.dk>
Date: Tue, 15 Feb 2022 07:36:38 +0000
Subject: [PATCH] add dAEP_dxy as cost_gradient_function in
 PyWakeAEPCostModelComponent

---
 examples/docs/example_2_wake_comparison.py |  3 ++-
 topfarm/cost_models/py_wake_wrapper.py     | 30 +++++++++++++++++-----
 2 files changed, 25 insertions(+), 8 deletions(-)

diff --git a/examples/docs/example_2_wake_comparison.py b/examples/docs/example_2_wake_comparison.py
index 631f25d9..5de212f4 100644
--- a/examples/docs/example_2_wake_comparison.py
+++ b/examples/docs/example_2_wake_comparison.py
@@ -58,7 +58,8 @@ def main():
         def get_tf(windFarmModel):
             return TopFarmProblem(
                 design_vars=dict(zip('xy', init_pos.T)),
-                cost_comp=PyWakeAEPCostModelComponent(windFarmModel, n_wt=3, ws=10, wd=np.arange(0, 360, 12)),
+                cost_comp=PyWakeAEPCostModelComponent(windFarmModel, n_wt=3, ws=10, wd=np.arange(0, 360, 12),
+                                                      grad_method=None),
                 constraints=[SpacingConstraint(min_spacing),
                              XYBoundaryConstraint(boundary)],
                 driver=EasyScipyOptimizeDriver(),
diff --git a/topfarm/cost_models/py_wake_wrapper.py b/topfarm/cost_models/py_wake_wrapper.py
index 800dc0f3..9dab4842 100644
--- a/topfarm/cost_models/py_wake_wrapper.py
+++ b/topfarm/cost_models/py_wake_wrapper.py
@@ -13,6 +13,7 @@ from scipy.interpolate.interpolate import RegularGridInterpolator
 import warnings
 from py_wake.flow_map import HorizontalGrid
 import xarray as xr
+from py_wake.utils.gradients import autograd
 
 
 class PyWakeAEP(AEPCalculator):
@@ -55,17 +56,32 @@ class PyWakeAEP(AEPCalculator):
 
 
 class PyWakeAEPCostModelComponent(AEPCostModelComponent):
-    def __init__(self, windFarmModel, n_wt, wd=None, ws=None, max_eval=None, **kwargs):
+    def __init__(self, windFarmModel, n_wt, wd=None, ws=None, max_eval=None, grad_method=autograd, **kwargs):
         self.windFarmModel = windFarmModel
+
+        def aep(**kwargs):
+            return self.windFarmModel.aep(x=kwargs[topfarm.x_key],
+                                          y=kwargs[topfarm.y_key],
+                                          h=kwargs.get(topfarm.z_key, None),
+                                          type=kwargs.get(topfarm.type_key, 0),
+                                          wd=wd, ws=ws)
+
+        if grad_method:
+            dAEPdxy = self.windFarmModel.dAEPdxy(grad_method)
+
+            def daep(**kwargs):
+                return dAEPdxy(x=kwargs[topfarm.x_key],
+                               y=kwargs[topfarm.y_key],
+                               h=kwargs.get(topfarm.z_key, None),
+                               type=kwargs.get(topfarm.type_key, 0),
+                               wd=wd, ws=ws)
+        else:
+            daep = None
         AEPCostModelComponent.__init__(self,
                                        input_keys=[topfarm.x_key, topfarm.y_key],
                                        n_wt=n_wt,
-                                       cost_function=lambda **kwargs: self.windFarmModel.aep(
-                                           x=kwargs[topfarm.x_key],
-                                           y=kwargs[topfarm.y_key],
-                                           h=kwargs.get(topfarm.z_key, None),
-                                           type=kwargs.get(topfarm.type_key, 0),
-                                           wd=wd, ws=ws),
+                                       cost_function=aep,
+                                       cost_gradient_function=daep,
                                        output_unit='GWh',
                                        max_eval=max_eval, **kwargs)
 
-- 
GitLab