(feat) Added MDAO api for aerodynamics optimization

The api features 3 components: The solver, the geometry and the initial mesh. The api is tested in a aerodynamic shape case where the NACA0012 is constrained at 0deg AoA with cl=0.40 and the objective function to minimize is the drag wrt the max camber and the max camber position. The test case ends with a final shape close to the NACA4412 (camber = 35% & max camber pos = 41%)

blast/api/mdaAPI.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2024 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+## BLASTER api for openMDAO
+# Paul Dechamps
+
+import numpy as np
+import openmdao.api as om
+import blast
+import blast.utils as viscUtils
+
+# Solver openMDAO object
+class BlasterSolver(om.ExplicitComponent):
+    """OpenMDAO component for the Blaster solver
+
+    Parameters
+    ----------
+    - cfg {dict} -- Configuration dictionary for the VII coupling
+    - icfg {dict} -- Configuration dictionary for the inviscid solver interface
+    - iSolverName {str} -- Inviscid solver name (default: 'DART')
+
+    Inputs
+    ------
+    - aoa {float} -- Angle of attack
+    - x_aero {np.array} -- Mesh coordinates
+
+    Outputs
+    -------
+    - cl {float} -- Lift coefficient
+    - cd {float} -- Drag coefficient
+    """
+
+    def __init__(self, cfg, icfg, iSolverName='DART'):
+        super().__init__()
+
+        self.nCpt = 0
+
+        if iSolverName == 'DART':
+            from blast.interfaces.dart.DartInterface import DartInterface as interface
+        else:
+            raise RuntimeError('Solver '+iSolverName+' currently not implemented')
+
+        for i in range(len(cfg['xtrF'])):
+            if cfg['xtrF'][i] is None:
+                cfg['xtrF'][i] = -1
+        self.vsol = blast.Driver(cfg['Re'], cfg['Minf'], cfg['CFL0'], cfg['nSections'], xtrF_top=cfg['xtrF'][0], xtrF_bot=cfg['xtrF'][1], _span=cfg['span'], _verbose=cfg['Verb'])
+
+        self.isol = interface(icfg, self.vsol, cfg)
+
+        # Coupler
+        import blast.coupler as blastCoupler
+        self.coupler = blastCoupler.Coupler(self.isol, self.vsol, _maxCouplIter=cfg['couplIter'], _couplTol=cfg['couplTol'], _iterPrint=cfg['iterPrint'], _resetInv=cfg['resetInv'])
+
+        # Adjoint
+        self.adjointSolver = blast.CoupledAdjoint(self.isol.adjointSolver, self.vsol)
+
+        # Single coupler run to initialize viscous sections (which depend on size of upper/lower sides)
+        self.coupler.run(write=False)
+        self.coupler.reset()
+        print('object created')
+    
+    def setup(self):
+        self.add_input('aoa', val=0.0)
+        self.add_input('x_aero', val=np.zeros(self.isol.iobj['bnd'].nodes.size() * self.isol.getnDim()))
+        self.add_output('cl', val=0.0)
+        self.add_output('cd', val=0.0)
+    
+    def setup_partials(self):
+        self.declare_partials('cl', 'aoa', method='exact')
+        self.declare_partials('cd', 'aoa', method='exact')
+        self.declare_partials('cl', 'x_aero', method='exact')
+        self.declare_partials('cd', 'x_aero', method='exact')
+
+
+    def compute_partials(self, inputs, partials, discrete_inputs=None):
+        """Compute partial derivatives using the adjoint solver
+        """
+        # Run adjoint solver
+        self.isol.iobj['adj'].run()
+        self.adjointSolver.reset()
+        self.adjointSolver.run()
+
+        # AoA derivatives
+        partials['cl', 'aoa'] = self.adjointSolver.tdCl_AoA
+        partials['cd', 'aoa'] = self.adjointSolver.tdCd_AoA
+
+        # Mesh derivative
+        nDim = self.isol.getnDim()
+        # Get mesh matrix
+        for inod, n in enumerate(self.isol.iobj['bnd'].nodes):
+            for idim in range(nDim):
+                partials['cl', 'x_aero'][0][inod*nDim + idim] = self.adjointSolver.tdCl_x[n.row][idim]
+                partials['cd', 'x_aero'][0][inod*nDim + idim] = self.adjointSolver.tdCd_x[n.row][idim]
+    
+    def compute(self, inputs, outputs):
+        """Update the mesh, aoa and run the solver
+        """
+        nDim = self.isol.getnDim()
+        # Modify surface and deform mesh
+        self.isol.iobj['mrf'].savePos()
+        for inod, n in enumerate(self.isol.iobj['bnd'].nodes):
+            for idim in range(nDim):
+                n.pos[idim] = inputs['x_aero'][inod*nDim + idim]
+        self.isol.iobj['mrf'].deform()
+
+        # Update viscous mesh
+        self.isol.updateMesh()
+
+        # Update angle of attack
+        aoa = inputs['aoa'][0]
+        self.isol.setAoA(aoa*np.pi/180) # aoa is given in deg
+
+        # Reset coupler
+        self.coupler.reset()
+
+        # Run forward problem
+        self.coupler.run(write=False)
+        outputs['cl'] = self.isol.getCl()
+        outputs['cd'] = self.isol.getCd() + self.vsol.Cdf
+
+        # Write results
+        self.isol.writeCp(sfx='_'+str(self.nCpt))
+        viscUtils.getSolution(self.isol.sec, write=True, toW=['x', 'y', 'cf'], sfx='_'+str(self.nCpt))
+        self.nCpt += 1
+
+    def getBoundaryMesh(self):
+        """Initial surface mesh coordinates
+        """
+        return BlasterMesh(nDim=self.isol.getnDim(), bnd=self.isol.iobj['bnd'])
+    
+    def getBoundaryGeometry(self):
+        """Airfoil geometry
+        """
+        return BlasterGeometry(nDim=self.isol.getnDim(), nNodes=self.isol.iobj['bnd'].nodes.size(), nodeRows=self.isol.vBnd[0][0].connectListNodes)
+
+# Surface mesh
+class BlasterMesh(om.IndepVarComp):
+    """Initial 2D surface mesh coordinates
+    
+    Parameters
+    ----------
+    - nDim {int} -- Number of dimensions
+    - bnd {object} -- Boundary object containing the mesh
+
+    Outputs
+    -------
+    - x_aero0 {np.array} -- Initial mesh coordinates
+    """
+
+    def initialize(self):
+        self.options.declare('nDim')
+        self.options.declare('bnd', recordable=False)
+    
+    def setup(self):
+        bnd = self.options['bnd']
+        nDim = self.options['nDim']
+        x_aero0 = np.zeros(bnd.nodes.size() * nDim)
+        for inod, n in enumerate(bnd.nodes):
+            for idim in range(nDim):
+                x_aero0[inod*nDim + idim] = n.pos[idim]
+        self.add_output('x_aero0', val=x_aero0)
+
+class BlasterGeometry(om.ExplicitComponent):
+    """Airfoil parametrization as a NACA 4-digit airfoil.
+    Uses the NACA 4-digit airfoil parametrization to generate an airfoil from a given x
+    distribution.
+    Parameters become tau (thickness), eps (max camber) and p (max camber position).
+
+    Parameters
+    ----------
+    - nDim {int} -- Number of dimensions
+    - nNodes {int} -- Number of nodes
+    - nodeRows {np.array} -- Nodal indexes
+
+    Inputs
+    ------
+    - x_aero0 {np.array} -- Initial mesh
+    - tau {float} -- Airfoil thickness
+    - eps {float} -- Airfoil max camber
+    - p {float} -- Airfoil max camber position
+
+    Outputs
+    -------
+    - x_aero_out {np.array} -- Airfoil coordinates
+    """
+
+    def initialize(self):
+        self.options.declare('nDim')
+        self.options.declare('nNodes')
+        self.options.declare('nodeRows')
+
+    def setup(self):
+        nNodes = self.options['nNodes']
+        nDim = self.options['nDim']
+        self.add_input('x_aero0', np.zeros(nNodes*nDim)) # Initial mesh
+        self.add_input('tau', val = 0.12, desc='Airfoil thickness')
+        self.add_input('eps', val = 0., desc='Airfoil max camber')
+        self.add_input('p', val = 0., desc='Airfoil max camber position')
+        self.add_output('x_aero_out', val = np.zeros(nNodes*nDim))
+
+    def setup_partials(self):
+        self.declare_partials('x_aero_out', 'tau', method='fd')
+        self.declare_partials('x_aero_out', 'eps', method='fd')
+        self.declare_partials('x_aero_out', 'p', method='fd')
+    
+    def compute(self, inputs, outputs):
+        """Compute airfoil coordinates
+        Points are sorted in seilig format to identify upper and lower surfaces.
+        Then the NACA 4-digit airfoil is generated and the points are sorted back
+        to the original mesh format."""
+
+        x_aero0 = inputs['x_aero0']
+        tau = inputs['tau'][0]
+        eps = inputs['eps'][0]
+        p = inputs['p'][0]
+        nDim = self.options['nDim']
+
+        # Retreive selig coordinates
+        x_seilig = np.zeros(len(self.options['nodeRows']))
+        for i, val in enumerate(self.options['nodeRows']):
+            x_seilig[i] = x_aero0[val*nDim]
+        # Generate airfoil
+        # Format tau eps and p with 2 digits after coma
+        print(f"tau: {tau:.10f}, eps: {eps:.10f}, p: {p:.10f}")
+        xb, yb = self.__generateAirfoil(x_seilig, tau, eps, p)
+        # Go back to normal mesh
+        x_aero_out = np.zeros(len(x_aero0))
+        for i, val in enumerate(self.options['nodeRows']):
+            x_aero_out[val*nDim] = xb[i]
+            x_aero_out[val*nDim+1] = yb[i]
+        outputs['x_aero_out'] = x_aero_out
+    
+    def __generateAirfoil(self, x, tau, eps, p):
+        """Generate NACA 4-digit airfoil from x coordinates points
+        
+        Arguments
+        ---------
+        - x {np.array} -- x coordinates sorted in seilig format
+        - tau {float} -- Airfoil thickness
+        - eps {float} -- Airfoil max camber
+        - p {float} -- Airfoil max camber position
+        
+        Returns
+        -------
+        - tuple -- x and y coordinates
+        """
+
+        T = 10 * tau * (
+          0.2969 * np.sqrt(x) 
+        - 0.1260 * (x) 
+        - 0.3537 * (x)**2 
+        + 0.2843 * (x)**3 
+        - 0.1015 * (x)**4) # Thickness
+
+        YBar     = np.zeros(len(T))
+        dYBar_dx = np.zeros(len(T))
+        Y        = np.zeros(len(T))
+        xB       = np.zeros(len(T))
+        yB       = np.zeros(len(T))
+
+        for i, xi in enumerate(x):
+            # Camber
+            if xi < p:
+                YBar[i]     = ((eps * xi)/(p**2)) * (2 * p - xi)
+                dYBar_dx[i] = (2 * eps)/(p**2) * (p - xi)
+            elif xi >= p:
+                YBar[i]     = (eps * (1 - xi)/((1-p)**2)) * (1 + xi - 2 * p)
+                dYBar_dx[i] = (2 * eps)/((1 - p)**2) * (p - xi)
+
+            # Final x & y coordinates
+            if i >= len(x)/2 - 1: 
+                # Lower surface
+                yB[i] = YBar[i] - 0.5 * T[i] * np.cos(np.arctan(dYBar_dx[i]))
+                xB[i] = x[i]    + 0.5 * T[i] * np.sin(np.arctan(dYBar_dx[i]))
+            else: 
+                # Upper surface
+                yB[i] = YBar[i] + 0.5 * T[i] * np.cos(np.arctan(dYBar_dx[i]))
+                xB[i] = x[i]    - 0.5 * T[i] * np.sin(np.arctan(dYBar_dx[i]))    
+        return (xB, yB)

blast/mdao/naca0012_mdao.py

+import blast.utils as viscUtils
+import numpy as np
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+from blast.api.mdaAPI import BlasterSolver
+
+import openmdao.api as om
+
+import math
+
+def cfgInviscid(nthrds, verb):
+    import os.path
+    # Parameters
+    return {
+    # Options
+    'scenario' : 'aerodynamic',
+    'task' : 'optimization',
+    'Threads' : nthrds, # number of threads
+    'Verb' : verb, # verbosity
+    # Model (geometry or mesh)
+    'File' : os.path.dirname(os.path.dirname(os.path.abspath(__file__))) + '/models/dart/n0012.geo', # Input file containing the model
+    'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF' : 30, 'msTe' : 0.01, 'msLe' : 0.0075}, # parameters for input file model
+    'Dim' : 2, # problem dimension
+    'Format' : 'gmsh', # save format (vtk or gmsh)
+    # Markers
+    'Fluid' : 'field', # name of physical group containing the fluid
+    'Farfield' : ['upstream', 'farfield', 'downstream'], # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
+    'Wing' : 'airfoil', # LIST of names of physical groups containing the lifting surface boundary
+    'Wake' : 'wake', # LIST of names of physical group containing the wake
+    'WakeTip' : 'wakeTip', # LIST of names of physical group containing the edge of the wake
+    'Te' : 'te', # LIST of names of physical group containing the trailing edge
+    'dbc' : False,
+    'Upstream' : 'upstream',
+    # Freestream
+    'M_inf' : 0.2, # freestream Mach number
+    'AoA' : 2., # freestream angle of attack
+    # Geometry
+    'S_ref' : 1., # reference surface length
+    'c_ref' : 1., # reference chord length
+    'x_ref' : .25, # reference point for moment computation (x)
+    'y_ref' : 0.0, # reference point for moment computation (y)
+    'z_ref' : 0.0, # reference point for moment computation (z)
+    # Numerical
+    'LSolver' : 'GMRES', # inner solver (Pardiso, MUMPS or GMRES)
+    'G_fill' : 2, # fill-in factor for GMRES preconditioner
+    'G_tol' : 1e-5, # tolerance for GMRES
+    'G_restart' : 50, # restart for GMRES
+    'Rel_tol' : 1e-6, # relative tolerance on solver residual
+    'Abs_tol' : 1e-8, # absolute tolerance on solver residual
+    'Max_it' : 20 # solver maximum number of iterations
+    }
+
+def cfgBlast(verb):
+    return {
+        'Re' : 1e7,                 # Freestream Reynolds number
+        'Minf' : 0.2,               # Freestream Mach number (used for the computation of the time step only)
+        'CFL0' : 1.,                 # Inital CFL number of the calculation
+        'Verb': verb,               # Verbosity level of the solver
+        'couplIter': 50,            # Maximum number of iterations
+        'couplTol' : 1e-3,          # Tolerance of the VII methodology
+        'iterPrint': 5,             # int, number of iterations between outputs
+        'resetInv' : True,          # bool, flag to reset the inviscid calculation at every iteration.
+        'sections' : [0],           # List of sections for boundary layer calculation
+        'xtrF' : [0.05, 0.05],       # Forced transition location
+        'interpolator' : 'Matching',# Interpolator for the coupling
+        'nSections' : 1,            # Number of sections in the domain
+        'span' : 1.0                # Wing span
+    }
+
+class groupMDA(om.Group):
+    def setup(self):
+        args = parseargs()
+        icfg = cfgInviscid(args.k, args.verb)
+        vcfg = cfgBlast(args.verb)
+
+        #cycle = self.add_subsystem('cycle', om.Group())
+
+        # Add the MDAO components
+        blaster = BlasterSolver(vcfg, icfg)
+
+        # Design variables box
+        self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*'])
+        self.add_subsystem('mesh', blaster.getBoundaryMesh())
+        self.add_subsystem('geometry', blaster.getBoundaryGeometry())
+        self.add_subsystem('solver', blaster)
+
+        # Add objective function
+        self.add_subsystem('obj_cmp', om.ExecComp('obj = cd'), promotes=['obj', 'cd'])
+
+        # Connect meshes
+        self.connect('mesh.x_aero0', 'geometry.x_aero0')
+        self.connect('geometry.x_aero_out', 'solver.x_aero')
+
+        #nonlinear_solver = om.NewtonSolver(solve_subsystems=False)
+        self.linear_solver = om.DirectSolver()
+    
+    def configure(self):
+        # Design vars
+        self.dvs.add_output('aoa', val=0.0)
+        self.connect('aoa', 'solver.aoa')
+        #self.add_design_var('aoa', lower=0.0, upper=3.0)
+
+        self.dvs.add_output('tau', val=0.12)
+        self.connect('tau', 'geometry.tau')
+        #self.add_design_var('tau', lower=0.10, upper=0.2)
+
+        self.dvs.add_output('eps', val=0.0)
+        self.connect('eps', 'geometry.eps')
+        self.add_design_var('eps', lower=0.0, upper=0.05)
+
+        self.dvs.add_output('p', val=0.4)
+        self.connect('p', 'geometry.p')
+        self.add_design_var('p', lower=0.35, upper=0.45)
+
+        self.add_constraint('solver.cl', equals=0.4)
+
+        #self.connect('solver.cl', 'cl')
+        self.connect('solver.cd', 'cd')
+        self.add_objective('obj')
+
+def main():
+    # Timer.
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    prob = om.Problem()
+    prob.model = groupMDA()
+
+    prob.driver = om.ScipyOptimizeDriver()
+    prob.driver.options['optimizer'] = 'SLSQP'
+    prob.driver.options['tol'] = 1e-3
+
+    recorder = om.SqliteRecorder('reports/BlasterSolver_file.sql')
+
+    prob.driver.add_recorder(recorder)
+    prob.driver.recording_options['record_objectives'] = True
+    prob.driver.recording_options['record_desvars'] = True
+    prob.driver.recording_options['includes'] = ['aoa', 'solver.cl', 'solver.cd', 'mesh.x_aero0', 'geometry.x_aero_out']
+
+    prob.setup()
+    #prob.check_partials(includes=['cl', 'aoa'], compact_print=True)
+    prob.run_driver()
+
+    # Instantiate your CaseReader
+    cr = om.CaseReader(prob.get_reports_dir() / "../BlasterSolver_file.sql")
+    # Isolate "problem" as your source
+    driver_cases = cr.list_cases('driver', out_stream=None)
+
+    # Get evolution of f and x
+    aoa_evol = []
+    cl_evol = []
+    cd_evol = []
+    shape = []
+    shape0 = []
+
+    for c in driver_cases:
+        case = cr.get_case(c)
+        aoa_evol.append(case.get_val('aoa')[0]*180/np.pi)
+        cl_evol.append(case.get_val('solver.cl')[0])
+        cd_evol.append(case.get_val('solver.cd')[0])
+        shape0.append(case.get_val('mesh.x_aero0'))
+        shape.append(case.get_val('geometry.x_aero_out'))
+
+    print('success')
+    print('aoa:', prob['aoa']*np.pi/180)
+    print('cl:', prob['solver.cl'])
+    print('cd:', prob['solver.cd'])
+    print('obj:', prob['obj'])
+
+    # Recover selig format
+    shape0_selig = np.array(shape0[0])
+    shape1_selig = np.array(shape[-1])
+    for i, val in enumerate(prob.model.solver.isol.vBnd[0][0].connectListNodes):
+        for j in range(prob.model.solver.isol.getnDim()):
+            shape0_selig[i*prob.model.solver.isol.getnDim()+j] = shape0[0][val*prob.model.solver.isol.getnDim()+j]
+            shape1_selig[i*prob.model.solver.isol.getnDim()+j] = shape[-1][val*prob.model.solver.isol.getnDim()+j]
+    from matplotlib import pyplot as plt
+
+    plt.plot(shape0_selig[0::2], shape0_selig[1::2], '--', color='grey', label='Initial shape - cl = {:.2f} - cd = {:.5f}'.format(cl_evol[0], cd_evol[0]))
+    plt.plot(shape1_selig[0::2], shape1_selig[1::2], '-', color='blue', label='Final shape - cl = {:.2f} - cd = {:.5f}'.format(cl_evol[-1], cd_evol[-1]))
+    plt.legend()
+    plt.axis('equal')
+    plt.show()
+    quit()
+
+
+
+    plt.plot(shape0[0][::2], shape0[0][1::2], 'x', color='grey')
+    for xx in shape:
+        plt.plot(xx[::2], xx[1::2], 'x', color='blue')
+    plt.axis('equal')
+    plt.show()
+    quit()
+
+
+if __name__ == "__main__":
+    main()