Merge pull request #26 from johnjasa/add_snopt

Add basic OpenMDAO and SNOPT wrapper for optimization

.github/workflows/CI_every_PR.yml

@@ -32,8 +32,20 @@ jobs:
           bash ./install.sh test-environment
           conda activate test-environment
         
-      - name: Run demo tests within WindSE
+      - name: Run all tests within WindSE
         shell: pwsh
         run: |
           conda activate test-environment
-          pytest -v --cov=windse tests/test_demos.py
+          pytest -v --cov=windse tests
+      
+      # Run coveralls
+      - name: Run coveralls
+        if: contains( matrix.os, 'ubuntu')
+        # This also works, https://github.com/AndreMiras/coveralls-python-action
+        #uses: AndreMiras/coveralls-python-action@develop
+        shell: pwsh
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+        run: |
+          conda activate test-environment
+          coveralls --service=github

.github/workflows/CI_every_commit.yml

@@ -1,7 +1,7 @@
 name: CI_every_commit
 
 # We run CI on push commits on all branches
-on: [push, pull_request]
+on: [push]
 
 # A workflow run is made up of one or more jobs that can run sequentially or in parallel
 jobs:
@@ -32,13 +32,6 @@ jobs:
           bash ./install.sh test-environment
           conda activate test-environment
 
-      # List the collected tests for debugging purposes
-      - name: List tests
-        shell: pwsh
-        run: |
-          conda activate test-environment
-          pytest --co
-        
       # Run all tests within WindSE, but not computationally expensive examples
       - name: Run tests within WindSE
         shell: pwsh

demo/documented/Yaml_Examples/7-wind_farm_2D_OM.yaml

@@ -0,0 +1,109 @@
+# General options
+general: 
+    name: 2_5D_Wind_Farm_Layout      # Name of the output folder
+    output: ["mesh","initial_guess","height","turbine_force","solution"]
+    dolfin_adjoint: true
+
+# Wind Farm Parameters: Uncomment a set to change the type of wind farm
+wind_farm: 
+
+    # ####################### Imported Wind Farm #######################
+    # type: imported                   #                          |
+    # path: Input_Data/wind_farm.txt   # location of wind farm    | -
+    # force: constant                      # radial force distribution| -
+    # ###################################################################
+
+
+    ######################### Grid Wind Farm #########################
+    type: grid            #                          |
+    grid_rows: 1          # Number of rows           | -
+    grid_cols: 2          # Number of columns        | -
+    ex_x: [-882, 882]     # x-extent of the farm     | m
+    ex_y: [-882, 882]     # y-extent of the farm     | m
+    HH: 90                # Hub Height               | m
+    RD: 126.0             # Turbine Diameter         | m
+    thickness: 12.0       # Effective Thickness      | m
+    yaw: 0.0              # Yaw                      | rads
+    axial: 0.25           # Axial Induction          | -
+    jitter: 0            # Randomly perturb turbines| m
+    seed: 8675309         # random seed for repeats  | -
+    force: constant       # radial force distribution| -
+    turbine_method: dolfin
+
+
+    ##################################################################
+
+
+    ########################## Grid Wind Farm #########################
+    # type: random          #                          |
+    # numturbs: 9           # number of turbines       | -
+    # ex_x: [-600, 600]     # x-extent of the farm     | m
+    # ex_y: [-600, 600]     # y-extent of the farm     | m
+    # seed: 5555555         # random seed for repeats  | -
+    # HH: 90                # Hub Height               | m
+    # RD: 126.0             # Turbine Diameter         | m
+    # thickness: 10.5       # Effective Thickness      | m
+    # yaw: 0.0              # Yaw                      | rads
+    # axial: 0.33           # Axial Induction          | -
+    # force: sine           # radial force distribution| -
+    ###################################################################
+
+
+
+# Domain Parameters: Uncomment a set to change domain shape
+domain: 
+
+    # ####################### Rectangle Domain #########################
+    # type: rectangle                              |
+    # x_range: [-1200, 1200] # x-range of the domain | m
+    # y_range: [-400, 400] # y-range of the domain | m
+    # nx: 300                # Number of x-nodes     | -
+    # ny: 100                # Number of y-nodes     | -
+    # ##################################################################
+
+
+    ########################### Circle Domain #########################
+    type: circle                                  
+    mesh_type: mshr         # squircular, elliptic, stretch
+    radius: 1500            # x-range of the domain | m
+    center: [0.0, 0.0]      # y-range of the domain | m
+    nt: 20                  # segments around circle| -
+    res: 10                 # resolution for mshr   | -
+    ###################################################################
+
+
+
+refine:
+    #                      # Description                  | Units
+    # farm_num: 1            # number of farm refinements   | -
+    # farm_type: square      # type of refinement at farm   | -
+    # farm_factor: 1.25       # farm radius multiplier       | -
+    turbine_num: 1         # number of turbine refinements| -
+    turbine_factor: 1.25   # turbine radius multiplier    | -
+
+function_space:
+    type: taylor_hood
+    turbine_degree: 6
+    turbine_space: Quadrature
+
+boundary_conditions:
+    vel_profile: uniform
+    HH_vel: 8.0
+
+problem:
+    use_25d_model: True
+    type: taylor_hood
+    viscosity: 5
+    lmax: 50
+
+solver:
+    nonlinear_solver: newton
+    newton_relaxation: 0.9
+    type: steady
+    save_power: true
+
+optimization:
+    control_types: [layout]
+    layout_bounds: [[-720, 720],[-720, 720]] 
+    optimize: True
+    opt_routine: OM_SLSQP

install.sh

@@ -14,7 +14,8 @@ conda install -y -c conda-forge fenics=2019.1.0=py38_9 dolfin-adjoint matplotlib
 pip install git+https://github.com/blechta/[email protected]
 pip install git+https://github.com/blechta/[email protected]
 pip install git+https://github.com/blechta/[email protected]
-pip install singledispatch networkx pulp
+pip install git+https://github.com/mdolab/[email protected]
+pip install singledispatch networkx pulp openmdao
 
 ### Install editible version of WindSE
 pip install -e .

windse/OptimizationManager.py

@@ -44,7 +44,141 @@
     import matplotlib.pyplot as plt
 else:
     InequalityConstraint = object
+
+
+import openmdao.api as om
+
+
+class ObjComp(om.ExplicitComponent):
+    def initialize(self):
+        self.options.declare('initial_DVs', types=np.ndarray)
+        self.options.declare('J', types=object)
+        self.options.declare('dJ', types=object)
+        self.options.declare('callback', types=object)
+
+    def setup(self):
+        self.add_input('DVs', val=self.options['initial_DVs'])
+        self.add_output('obj', val=0.)
+
+        self.declare_partials('*', '*')
+
+    def compute(self, inputs, outputs):
+        m = list(inputs['DVs'])
+        computed_output = self.options['J'](m)
+        outputs['obj'] = computed_output
+        self.options['callback'](m)
+
+    def compute_partials(self, inputs, partials):
+        m = list(inputs['DVs'])
+        jac = self.options['dJ'](m)
+        partials['obj', 'DVs'] = jac
+
+
+class ConsComp(om.ExplicitComponent):
+    def initialize(self):
+        self.options.declare('initial_DVs', types=np.ndarray)
+        self.options.declare('J', types=object)
+        self.options.declare('dJ', types=object)
+        self.options.declare('con_name', types=str)
+
+    def setup(self):
+        self.con_name = self.options["con_name"]
+        self.add_input('DVs', val=self.options['initial_DVs'])
+
+        output = self.options['J'](self.options['initial_DVs'])
+        self.add_output(self.con_name, val=output)
 
+        self.declare_partials('*', '*')
+
+    def compute(self, inputs, outputs):
+        m = list(inputs['DVs'])
+        computed_output = self.options['J'](m)
+        outputs[self.con_name] = computed_output
+
+    def compute_partials(self, inputs, partials):
+        m = list(inputs['DVs'])
+        jac = self.options['dJ'](m)
+        partials[self.con_name, 'DVs'] = jac
+
+
+def gather(m):
+    if isinstance(m, list):
+        return list(map(gather, m))
+    elif hasattr(m, "_ad_to_list"):
+        return m._ad_to_list(m)
+    else:
+        return m  # Assume it is gathered already
+
+def om_wrapper(J, initial_DVs, dJ, H, bounds, **kwargs):
+
+    # build the model
+    prob = om.Problem(model=om.Group())
+
+    if 'callback' in kwargs:
+        callback = kwargs['callback']
+    else:
+        callback = None
+
+    prob.model.add_subsystem('obj_comp', ObjComp(initial_DVs=initial_DVs, J=J, dJ=dJ, callback=callback), promotes=['*'])
+
+    constraint_types = []
+    if 'constraints' in kwargs:
+        constraints = kwargs['constraints']
+
+        if not isinstance(constraints, list):
+            constraints = [constraints]
+
+        for idx, c in enumerate(constraints):
+            if isinstance(c, InequalityConstraint):
+                typestr = "ineq"
+            elif isinstance(c, EqualityConstraint):
+                typestr = "eq"
+            else:
+                raise Exception("Unknown constraint class")
+
+            def jac(x):
+                out = c.jacobian(x)
+                return [gather(y) for y in out]
+
+            constraint_types.append(typestr)
+
+            con_name = f'con_{idx}'
+            prob.model.add_subsystem(f'cons_comp_{idx}', ConsComp(initial_DVs=initial_DVs, J=c.function, dJ=jac, con_name=con_name), promotes=['*'])
+
+    lower_bounds = bounds[:, 0]
+    upper_bounds = bounds[:, 1]
+
+    # set up the optimization
+    if 'SLSQP' in kwargs['opt_routine']:
+        prob.driver = om.ScipyOptimizeDriver()
+        prob.driver.options['optimizer'] = 'SLSQP'
+    elif 'SNOPT' in kwargs['opt_routine']:
+        prob.driver = om.pyOptSparseDriver()
+        prob.driver.options['optimizer'] = 'SNOPT'
+
+    prob.model.add_design_var('DVs', lower=lower_bounds, upper=upper_bounds)
+    prob.model.add_objective('obj')
+
+    for idx, constraint_type in enumerate(constraint_types):
+        con_name = f'con_{idx}'
+        if constraint_type == "eq":
+            prob.model.add_constraint(con_name, equals=0.)
+        else:
+            # Inequality means it's positive from scipy and dolfin
+            prob.model.add_constraint(con_name, lower=0.)            
+
+    prob.setup()
+
+    prob.set_val('DVs', initial_DVs)
+
+    # Run the optimization
+    prob.run_driver()
+
+    # Return the optimal design variables
+    return(prob['DVs'])
+
+
+
 class Optimizer(object):
     """
     A GenericProblem contains on the basic functions required by all problem objects.
@@ -72,7 +206,7 @@ def __init__(self, solver):
 
         ### Process parameters ###
         if "layout" in self.control_types:
-            if isinstance(self.layout_bounds,list):
+            if isinstance(self.layout_bounds,(list, np.ndarray)):
                 self.layout_bounds = np.array(self.layout_bounds)
             elif self.layout_bounds == "wind_farm":
                 self.layout_bounds = np.array([self.farm.ex_x,self.farm.ex_y])
@@ -401,14 +535,20 @@ def get_minimum_distance_constraint_func(self, m_pos, min_distance=200):
     def Optimize(self):
 
         self.fprint("Beginning Optimization",special="header")
-
-        if "layout" in self.control_types:
-            self.m_opt=minimize(self.Jhat, method="SLSQP", options = {"disp": True}, constraints = self.dist_constraint, bounds = self.bounds, callback = self.OptPrintFunction)
+
+        # TODO : simplify this logic
+        if "SNOPT" in self.opt_routine or "OM_SLSQP" in self.opt_routine:
+            if "layout" in self.control_types:
+                m_opt=minimize(self.Jhat, method="Custom", options = {"disp": True}, constraints = self.dist_constraint, bounds = self.bounds, callback = self.OptPrintFunction, algorithm=om_wrapper, opt_routine=self.opt_routine)
+            else:
+                m_opt=minimize(self.Jhat, method="Custom", options = {"disp": True}, bounds = self.bounds, callback = self.OptPrintFunction, algorithm=om_wrapper, opt_routine=self.opt_routine)
         else:
-            # m_opt=minimize(self.Jhat, method="SLSQP", options = {"disp": True}, bounds = self.bounds,  callback = self.OptPrintFunction)
-            # m_opt=minimize(self.Jhat, method="L-BFGS-B", options = {"disp": True}, bounds = self.bounds, callback = self.OptPrintFunction)
-            self.m_opt=minimize(self.Jhat, method=self.opt_routine, options = {"disp": True}, bounds = self.bounds, callback = self.OptPrintFunction)
-
+            if "layout" in self.control_types:
+                m_opt=minimize(self.Jhat, method=self.opt_routine, options = {"disp": True}, constraints = self.dist_constraint, bounds = self.bounds, callback = self.OptPrintFunction)
+            else:
+                m_opt=minimize(self.Jhat, method=self.opt_routine, options = {"disp": True}, bounds = self.bounds, callback = self.OptPrintFunction)
+
+        self.m_opt = m_opt
 
         if self.num_controls == 1:
             self.m_opt = (self.m_opt,)