(feat) Advanced opti cases

+ Refactor adjoint and api's

blast/api/mdaAPI.py

@@ -20,9 +20,13 @@
 # Paul Dechamps
 
 import numpy as np
+from fwk.coloring import ccolors
 import openmdao.api as om
 import blast
 import blast.utils as viscUtils
+from pygeo import DVGeometryCST
+from scipy.optimize import least_squares
+import fwk
 
 # Solver openMDAO object
 class BlasterSolver(om.ExplicitComponent):
@@ -49,6 +53,7 @@ class BlasterSolver(om.ExplicitComponent):
         super().__init__()
 
         self.nCpt = 0
+        self.tms = fwk.Timers()
 
         if iSolverName == 'DART':
             from blast.interfaces.dart.DartInterface import DartInterface as interface
@@ -75,10 +80,11 @@ class BlasterSolver(om.ExplicitComponent):
         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)
+        self.add_input('aoa', val=0.0, desc='Flow angle of attack')
+        self.add_input('x_aero', val=np.zeros(self.isol.iobj['bnd'].nodes.size() * self.isol.getnDim()), desc='Surface mesh coordinates')
+        self.add_input('visc_mode', val=True) # Inviscid or viscous mode
+        self.add_output('cl', val=0.0, desc='Lift coefficient')
+        self.add_output('cd', val=0.0, desc='Drag coefficient')
     
     def setup_partials(self):
         self.declare_partials('cl', 'aoa', method='exact')
@@ -90,63 +96,177 @@ class BlasterSolver(om.ExplicitComponent):
     def compute_partials(self, inputs, partials, discrete_inputs=None):
         """Compute partial derivatives using the adjoint solver
         """
-        # Run adjoint solver
+        nDim = self.isol.getnDim()
+        self.tms['adj'].start()
+
+        print(ccolors.ANSI_BLUE, 'Computing derivatives in mode '+('viscous' if inputs['visc_mode'][0] else 'inviscid'), ccolors.ANSI_RESET)
+        
+        # Run inviscid adjoint solver
         self.isol.iobj['adj'].run()
-        self.adjointSolver.reset()
-        self.adjointSolver.run()
+
+        dClaoa = 0.
+        dCdaoa = 0.
+        dClx = np.zeros(self.isol.iobj['bnd'].nodes.size() * nDim)
+        dCdx = np.zeros(self.isol.iobj['bnd'].nodes.size() * nDim)
+
+        if inputs['visc_mode'][0]:
+            self.adjointSolver.reset()
+            self.adjointSolver.run()
+            dClaoa = self.adjointSolver.tdCl_AoA
+            dCdaoa = self.adjointSolver.tdCd_AoA
+
+            # Mesh derivatives
+            for inod, n in enumerate(self.isol.iobj['bnd'].nodes):
+                for idim in range(nDim):
+                    dClx[inod*nDim + idim] = self.adjointSolver.tdCl_x[n.row][idim]
+                    dCdx[inod*nDim + idim] = self.adjointSolver.tdCd_x[n.row][idim]
+        elif not inputs['visc_mode'][0]:
+            dClaoa = self.isol.iobj['adj'].dClAoa
+            dCdaoa = self.isol.iobj['adj'].dCdAoa
+
+            # Mesh derivatives
+            for inod, n in enumerate(self.isol.iobj['bnd'].nodes):
+                for idim in range(nDim):
+                    dClx[inod*nDim + idim] = self.isol.iobj['adj'].dClMsh[n.row][idim]
+                    dCdx[inod*nDim + idim] = self.isol.iobj['adj'].dCdMsh[n.row][idim]
+        else:
+            raise ValueError('omBlasterSolver::compute_partials - Invalid mode', inputs['visc_mode'][0])
 
         # AoA derivatives
-        partials['cl', 'aoa'] = self.adjointSolver.tdCl_AoA
-        partials['cd', 'aoa'] = self.adjointSolver.tdCd_AoA
+        partials['cl', 'aoa'] = dClaoa * np.pi/180 # /°
+        partials['cd', 'aoa'] = dCdaoa * np.pi/180 # /°
+        partials['cl', 'x_aero'][0] = dClx
+        partials['cd', 'x_aero'][0] = dCdx
+        self.tms['adj'].stop()
 
-        # 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):
+    def compute(self, inputs, outputs, discrete_inputs=None):
         """Update the mesh, aoa and run the solver
         """
+        self.tms['update'].start()
         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()
+        self.tms['update'].stop()
 
         # Run forward problem
-        self.coupler.run(write=False)
-        outputs['cl'] = self.isol.getCl()
-        outputs['cd'] = self.isol.getCd() + self.vsol.Cdf
+        self.tms['fwd'].start()
+        print(ccolors.ANSI_BLUE, 'Computing flow for AoA=', round(aoa,2), 'in mode '+('viscous' if inputs['visc_mode'][0] else 'inviscid'), ccolors.ANSI_RESET)
+
+        if inputs['visc_mode'][0]:
+            self.coupler.reset()
+            self.coupler.run(write=False)
+        elif not inputs['visc_mode'][0]:
+            self.isol.reset()
+            self.isol.run()
+        else:
+            raise ValueError('omBlasterSolver::compute - Invalid mode', inputs['visc_mode'][0])
+        self.tms['fwd'].stop()
+
+        # write solution
+        self.__write()
 
-        # Write results
-        self.isol.writeCp(sfx='_'+str(self.nCpt))
-        viscUtils.getSolution(self.isol.sec, write=True, toW=['x', 'y', 'cf'], sfx='_'+str(self.nCpt))
+        outputs['cl'] = self.isol.getCl()
+        outputs['cd'] = self.isol.getCd()
+        if inputs['visc_mode'][0]:
+            outputs['cd'] += self.vsol.Cdf # Add viscous contribution
         self.nCpt += 1
 
+    def __write(self):
+        saveAll = False
+        if self.nCpt == 0:
+            sfx = '_0'
+        elif saveAll:
+            sfx = '_'+str(self.nCpt)
+        else:
+            sfx = '_opt'
+        viscUtils.getSolution(self.isol.sec, write=True, toW=['x', 'y', 'cf'], sfx=sfx) # Boundary layer
+        self.isol.writeCp(sfx=sfx) # Cp
+        self.isol.iobj['wrt'].save(self.isol.iobj['pbl'].msh.name + sfx) # Mesh
+        self.isol.save(sfx) # Solution
+        np.savetxt('shape'+sfx+'.txt', self.isol.vBnd[0][0].nodesCoord) # Shape
+        np.savetxt('loads'+sfx+'.txt', [self.isol.getCl(), self.isol.getCd() + self.vsol.Cdf, self.isol.getCm()]) # Loads
+
     def getBoundaryMesh(self):
         """Initial surface mesh coordinates
         """
         return BlasterMesh(nDim=self.isol.getnDim(), bnd=self.isol.iobj['bnd'])
     
-    def getBoundaryGeometry(self):
-        """Airfoil geometry
+    def getBoundaryNACAGeometry(self):
+        """Airfoil naca-4digits geometry
+        """
+        return BlasterNACAGeometry(nDim=self.isol.getnDim(), nNodes=self.isol.iobj['bnd'].nodes.size(), nodeRows=self.isol.vBnd[0][0].connectListNodes)
+    
+    def getBoundaryFFDGeometry(self):
+        """Airfoil FFD geometry
         """
-        return BlasterGeometry(nDim=self.isol.getnDim(), nNodes=self.isol.iobj['bnd'].nodes.size(), nodeRows=self.isol.vBnd[0][0].connectListNodes)
+        with open('surface_DVGeo.dat', 'w') as f:
+            for i, node_coord in enumerate(self.isol.vBnd[0][0].nodesCoord):
+                f.write(f"{node_coord[0]} {node_coord[1]}\n")
+        
+        _dvGeo = DVGeometryCST('surface_DVGeo.dat', numCST=4)
+        _dvGeo.addDV("upper_shape", dvType="upper", lowerBound=-0.1, upperBound=0.5)
+        _dvGeo.addDV("lower_shape", dvType="lower", lowerBound=-0.5, upperBound=0.1)
+        _dvGeo.addDV("class_shape_n1", dvType="N1")
+        _dvGeo.addDV("class_shape_n2", dvType="N2")
+        _dvGeo.addDV("chord", dvType="chord")
+        points = np.array(np.loadtxt('surface_DVGeo.dat'))
+        points = np.hstack((points, np.zeros((points.shape[0], 1))))  # add 0s for z coordinates (unused)
+        _dvGeo.addPointSet(points, "airfoil_coords")
+        _dvGeo.update("airfoil_coords")
+        return BlasterFFDGeometry(nDim=self.isol.getnDim(), nNodes=self.isol.iobj['bnd'].nodes.size(), nodeRows=self.isol.vBnd[0][0].connectListNodes, dvGeo=_dvGeo)
+
+# Updater to activate viscous mode
+class BlasterModeUpdater(om.ExplicitComponent):
+    def initialize(self):
+        self.options.declare('tolerance', default=1e-6, types=float, desc='Tolerance value to trigger visc_mode')
+        self.options.declare('force_mode', default=None, types=str, desc='Force ''visc'' mode or ''inviscid'' mode')
+
+    def setup(self):
+        self.add_input('current_objective', val=1.0)
+        self.add_output('visc_mode', val=0.0)
+        self._previous_objective = 0.0
+        self._lock = False
+
+    def compute(self, inputs, outputs):
+
+        # Forced mode case
+        if self.options['force_mode'] is not None:
+            print('coucou', self.options['force_mode'])
+            if self.options['force_mode'] == 'visc':
+                print('passe dans visc')
+                outputs['visc_mode'] = True
+            elif self.options['force_mode'] == 'inv':
+                print('passe dans inv')
+                outputs['visc_mode'] = False
+            else:
+                raise RuntimeError('Invalid force_mode', self.options['force_mode'])
+            print('output visc1', outputs['visc_mode'])
+            return
+
+        tolerance = self.options['tolerance']
+        relDiff = abs((inputs['current_objective']-self._previous_objective)/inputs['current_objective'])
+        print('curr', inputs['current_objective'][0], 'prev', self._previous_objective)
+        print('relDiff', relDiff, 'tolerance', tolerance)
+        if relDiff <= tolerance and relDiff != 0. or self._lock:
+            outputs['visc_mode'] = True
+            self._lock = True
+        else:
+            outputs['visc_mode'] = False
+        print('output visc', outputs['visc_mode'])
+        self._previous_objective = inputs['current_objective'][0]
+    
+    def compute_partials(self, inputs, partials):
+        pass
 
 # Surface mesh
 class BlasterMesh(om.IndepVarComp):
@@ -175,7 +295,7 @@ class BlasterMesh(om.IndepVarComp):
                 x_aero0[inod*nDim + idim] = n.pos[idim]
         self.add_output('x_aero0', val=x_aero0)
 
-class BlasterGeometry(om.ExplicitComponent):
+class BlasterNACAGeometry(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.
@@ -292,3 +412,103 @@ class BlasterGeometry(om.ExplicitComponent):
                 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)
+
+class BlasterFFDGeometry(om.ExplicitComponent):
+    """Airfoil parametrization as a FFD airfoil.
+    Uses the FFD airfoil parametrization to generate an airfoil from a given x
+    distribution.
+    """
+
+    def initialize(self):
+        self.options.declare('nDim', desc='Number of dimensions')
+        self.options.declare('nNodes', desc='Number of surface nodes')
+        self.options.declare('nodeRows', desc='Nodal index of the surface nodes sorted in seilig format')
+        self.options.declare('dvGeo', desc='Parametrization of the geometry', recordable=False)
+
+    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('upper_shape', [0.17553542, 0.14723497, 0.14531762, 0.13826955])
+        self.add_input('lower_shape', [-0.16935506, -0.15405552, -0.13912672, -0.14168947])
+        self.add_input('class_shape_n1', 0.5)
+        self.add_input('class_shape_n2', 1.0)
+        self.add_output('x_aero_out', val = np.zeros(nNodes*nDim))
+        self.add_output('thickness', val = .12)
+        self.add_output('LEradius', val = 0.02)
+
+    def setup_partials(self):
+        self.declare_partials('x_aero_out', 'upper_shape', method='fd')
+        self.declare_partials('x_aero_out', 'lower_shape', method='fd')
+        #self.declare_partials('x_aero_out', 'class_shape_n1', method='fd')
+        #self.declare_partials('x_aero_out', 'class_shape_n2', method='fd')
+        self.declare_partials('thickness', 'upper_shape', method='fd')
+        self.declare_partials('thickness', 'lower_shape', method='fd')
+        # self.declare_partials('thickness', 'class_shape_n1', method='fd')
+        # self.declare_partials('thickness', 'class_shape_n2', method='fd')
+        self.declare_partials('LEradius', 'upper_shape', method='fd')
+        self.declare_partials('LEradius', 'lower_shape', method='fd')
+        # self.declare_partials('LEradius', 'class_shape_n1', method='fd')
+        # self.declare_partials('LEradius', 'class_shape_n2', method='fd')
+
+    def compute(self, inputs, outputs):
+        nDim = self.options['nDim']
+
+        self.options['dvGeo'].setDesignVars({"upper_shape": inputs['upper_shape'],
+                             "lower_shape": inputs['lower_shape'],
+                             "class_shape_n1": inputs['class_shape_n1'],
+                             "class_shape_n2": inputs['class_shape_n2']})
+        points1 = self.options['dvGeo'].update("airfoil_coords")
+
+        # Compute thickness at 1/4 chord
+        upper_side = points1[:len(points1)//2]
+        lower_side = points1[len(points1)//2:]
+        arg_up = np.argmin(np.abs(upper_side[:, 0] - 0.38))
+        arg_low = np.argmin(np.abs(lower_side[:, 0] - 0.38))
+        maxThickness = upper_side[arg_up, 1] - lower_side[arg_low, 1]
+
+        # Extract leading edge points (assuming points are sorted by x)
+        le_points = points1[np.abs(points1[:, 0]) < 0.01]  # Adjust threshold as needed
+
+        # Separate upper and lower surface points
+        upper_surface = le_points[le_points[:, 1] >= 0]
+        lower_surface = le_points[le_points[:, 1] < 0]
+
+        # Combine upper and lower surface points for circle fitting
+        le_points_combined = np.vstack((upper_surface, lower_surface))
+
+        # Fit circle to leading edge points
+        center, radius = self.fit_circle(le_points_combined[:, 0], le_points_combined[:, 1])
+
+        # print('Upper shape', inputs['upper_shape'])
+        # print('Lower shape', inputs['lower_shape'])
+        # print('Max thickness', round(maxThickness, 4))
+        # print('Center:', center)
+        # print('Radius:', radius)
+
+        # Go back to normal mesh
+        x_aero_out = np.zeros(len(inputs['x_aero0']))
+        for i, val in enumerate(self.options['nodeRows']):
+            for idim in range(nDim):
+                x_aero_out[val*nDim+idim] = points1[i, idim]
+        outputs['x_aero_out'] = x_aero_out
+        outputs['thickness'] = maxThickness
+        outputs['LEradius'] = radius
+
+    def fit_circle(self, x, y):
+        def calc_R(xc, yc):
+            return np.sqrt((x - xc)**2 + (y - yc)**2)
+
+        def f_2(c):
+            Ri = calc_R(*c)
+            return Ri - Ri.mean()
+
+        x_m = np.mean(x)
+        y_m = np.mean(y)
+        center_estimate = x_m, y_m
+        result = least_squares(f_2, center_estimate)
+        center = result.x
+        Ri = calc_R(*center)
+        R = Ri.mean()
+        return center, R

blast/mdao/naca0012_mdao.py → blast/mdao/opti_NACA.py

@@ -85,7 +85,7 @@ class groupMDA(om.Group):
         # 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('geometry', blaster.getBoundaryNACAGeometry())
         self.add_subsystem('solver', blaster)
 
         # Add objective function

blast/mdao/opti_subsonic.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' : 25, 'msTe' : 0.0075, 'msLe' : 0.001}, # 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.1, # freestream Mach number
+    'AoA' : 1., # 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' : 'SparseLu', # 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.1,               # 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-4,          # 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' : [None, None],       # 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.getBoundaryFFDGeometry())
+        # Add an intermediate variable for the constraint (such that upper_shape[0] = lower_shape[0])
+        self.add_subsystem('constraint_comp', om.ExecComp('constraint_var = lower_shape[0] + upper_shape[0]',
+                                                          lower_shape=np.zeros(4), upper_shape=np.zeros(4)), promotes=['*'])
+        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=1.0)
+        self.connect('aoa', 'solver.aoa')
+        #self.add_design_var('aoa', lower=0.0, upper=3.0)
+
+        self.dvs.add_output('upper_shape', val=[0.17553542, 0.14723497, 0.14531762, 0.13826955])
+        self.connect('upper_shape', 'geometry.upper_shape')
+        self.add_design_var('upper_shape', lower=-0.1, upper=0.4)
+
+        self.dvs.add_output('lower_shape', val=[-0.17553542, -0.15405552, -0.13912672, -0.14168947])
+        self.connect('lower_shape', 'geometry.lower_shape')
+        self.add_design_var('lower_shape', lower=-0.5, upper=0.2)
+
+        self.dvs.add_output('class_shape_n1', val=0.5)
+        self.connect('class_shape_n1', 'geometry.class_shape_n1')
+        # self.add_design_var('class_shape_n1', lower=-0.5, upper=0.1)
+
+        self.dvs.add_output('class_shape_n2', val=1.)
+        self.connect('class_shape_n2', 'geometry.class_shape_n2')
+        # self.add_design_var('class_shape_n2', lower=-0.5, upper=0.1)
+
+        # Constraints
+        # Cl constraint
+        self.add_constraint('solver.cl', equals=0.6)
+        self.add_constraint('geometry.thickness', equals=.09)
+        self.add_constraint('geometry.LEradius', equals=0.0224)
+        self.add_constraint('constraint_var', equals=0., linear=True)
+
+        #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()
+
+    print(prob.model.solver.tms)
+
+    # 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()

blast/mdao/opti_transonic.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
+from blast.api.mdaAPI import BlasterModeUpdater
+
+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' : 20, 'msTe' : 0.01, 'msLe' : 0.002}, # 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.1, # freestream Mach number
+    'AoA' : 0., # 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' : 'SparseLu', # 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' : 1e8,                 # Freestream Reynolds number
+        'Minf' : 0.1,               # 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': 75,            # Maximum number of iterations
+        'couplTol' : 1e-5,          # 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)
+        gr = 1.055
+        if gr == 1.0:
+            icfg['Pars']['msF'] = icfg['Pars']['msLe']
+        else:
+            n = np.log10(1-(1-gr)*icfg['Pars']['xLgt']/icfg['Pars']['msLe'])/np.log10(gr)
+            icfg['Pars']['msF'] = icfg['Pars']['msLe']*gr**(n-1)
+
+        # Add the MDAO components
+        blaster = BlasterSolver(vcfg, icfg)
+
+        self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*']) # Design variables box
+        self.add_subsystem('mod', BlasterModeUpdater(tolerance=10., force_mode='visc')) # Mode updater
+        self.add_subsystem('mesh', blaster.getBoundaryMesh()) # Initial mesh
+        self.add_subsystem('geometry', blaster.getBoundaryFFDGeometry()) # Geometry computation
+        # Add an intermediate variable for the constraint (such that upper_shape[0] = lower_shape[0])
+        self.add_subsystem('constraint_comp', om.ExecComp('constraint_var = lower_shape[0] + upper_shape[0]',
+                                                          lower_shape=np.zeros(4), upper_shape=np.zeros(4)), promotes=['*'])
+        self.add_subsystem('solver', blaster) # Blaster solver
+        self.add_subsystem('obj_cmp', om.ExecComp('obj = cd'), promotes=['obj', 'cd']) # Objective function
+
+        # Connect meshes
+        self.connect('mesh.x_aero0', 'geometry.x_aero0')
+        self.connect('geometry.x_aero_out', 'solver.x_aero')
+        # Connect mode updater
+        self.connect('mod.visc_mode', 'solver.visc_mode')
+
+        #nonlinear_solver = om.NewtonSolver(solve_subsystems=False)
+        self.linear_solver = om.DirectSolver()
+    
+    def configure(self):
+        # Design vars
+        self.dvs.add_output('aoa', val=0.)
+        self.connect('aoa', 'solver.aoa')
+        self.add_design_var('aoa', lower=-5.0, upper=5.0)
+
+        self.dvs.add_output('upper_shape', val=[0.17553542, 0.14723497, 0.14531762, 0.13826955]) # NACA 0012
+        #self.dvs.add_output('upper_shape', val=[0.25367349, 0.24936591, 0.26996167, 0.28365274]) # NACA4412
+        self.connect('upper_shape', 'geometry.upper_shape')
+        self.add_design_var('upper_shape', lower=-0.1, upper=0.22)
+
+        self.dvs.add_output('lower_shape', val=[-0.16935506, -0.15405552, -0.13912672, -0.14168947]) # NACA 0012
+        #self.dvs.add_output('lower_shape', val=[-0.1330126,  -0.01246078, -0.04065711, -0.0123221]) # NACA 4412
+        self.connect('lower_shape', 'geometry.lower_shape')
+        self.add_design_var('lower_shape', lower=-0.2, upper=0.2)
+
+        self.dvs.add_output('class_shape_n1', val=0.5)
+        self.connect('class_shape_n1', 'geometry.class_shape_n1')
+        # self.add_design_var('class_shape_n1', lower=-0.5, upper=0.1)
+
+        self.dvs.add_output('class_shape_n2', val=1.)
+        self.connect('class_shape_n2', 'geometry.class_shape_n2')
+        # self.add_design_var('class_shape_n2', lower=-0.5, upper=0.1)
+
+        #self.connect('obj', 'mod.current_objective')
+
+        # Constraints
+        # Cl constraint
+        self.add_constraint('solver.cl', equals=0.7)
+        self.add_constraint('geometry.thickness', lower=0.12)
+        self.add_constraint('geometry.LEradius', equals=0.0224)
+        self.add_constraint('constraint_var', equals=0., linear=True)
+
+        #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()
+    print('check')
+    #prob.check_partials(includes=['cl', 'cd'], compact_print=True)
+    print('done')
+    prob.run_driver()
+
+    print(prob.model.solver.tms)
+
+    # 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])
+        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'])
+    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()

blast/src/DCoupledAdjoint.cpp

@@ -349,8 +349,15 @@ void CoupledAdjoint::run()
     tms["2-Solve"].stop();
 
     // Print output
-    // aoa gradients
     std::cout << "-------------- Adjoint solution --------------" << std::endl;
+    // aoa gradients
+    std::cout << std::setw(15) << std::left << "Adjoint AoA"
+              << std::setw(15) << std::right << "Res[Cl_AoA]"
+              << std::setw(15) << std::right << "Res[Cd_AoA]" << std::endl;
+    std::cout << std::fixed << std::setprecision(5);
+    std::cout << std::setw(15) << std::left << ""
+              << std::setw(15) << std::right << log10((A_adjoint * Eigen::Map<Eigen::VectorXd>(lambdaCl.data(), lambdaCl.size()) - Eigen::Map<Eigen::VectorXd>(rhsCl.data(), rhsCl.size())).norm())
+              << std::setw(15) << std::right << log10((A_adjoint * Eigen::Map<Eigen::VectorXd>(lambdaCd.data(), lambdaCd.size()) - Eigen::Map<Eigen::VectorXd>(rhsCd.data(), rhsCd.size())).norm()) << std::endl;
     std::cout << std::setw(15) << std::left << "AoA gradients"
               << std::setw(15) << std::right << "dCl_dAoA"
               << std::setw(15) << std::right << "dCd_dAoA" << std::endl;
@@ -359,6 +366,13 @@ void CoupledAdjoint::run()
               << std::setw(15) << std::right << tdCl_AoA
               << std::setw(15) << std::right << tdCd_AoA << std::endl;
     // mesh gradients
+    std::cout << std::setw(15) << std::left << "Adjoint mesh"
+              << std::setw(15) << std::right << "Res[Cl_x]"
+              << std::setw(15) << std::right << "Res[Cd_x]" << std::endl;
+    std::cout << std::fixed << std::setprecision(5);
+    std::cout << std::setw(15) << std::left << ""
+              << std::setw(15) << std::right << log10((dRx_x.transpose() * Eigen::Map<Eigen::VectorXd>(lambdaCl_x.data(), lambdaCl_x.size()) - Eigen::Map<Eigen::VectorXd>(rhsCl_x.data(), rhsCl_x.size())).norm())
+              << std::setw(15) << std::right << log10((dRx_x.transpose() * Eigen::Map<Eigen::VectorXd>(lambdaCd_x.data(), lambdaCd_x.size()) - Eigen::Map<Eigen::VectorXd>(rhsCd_x.data(), rhsCd_x.size())).norm()) << std::endl;
     std::cout << std::setw(15) << std::left << "Mesh gradients"
               << std::setw(15) << std::right << "Norm[dCl_dX]"
               << std::setw(15) << std::right << "Norm[dCd_dX]" << std::endl;
@@ -366,7 +380,7 @@ void CoupledAdjoint::run()
     std::cout << std::setw(15) << std::left << ""
               << std::setw(15) << std::right << Eigen::Map<Eigen::VectorXd>(lambdaCl_x.data(), lambdaCl_x.size()).norm()
               << std::setw(15) << std::right << Eigen::Map<Eigen::VectorXd>(lambdaCd_x.data(), lambdaCd_x.size()).norm() << std::endl;
-
+    tms["0-Total"].stop();
     std::cout << "--------------- Adjoint timers ---------------" << std::endl;
     std::cout << tms << "----------------------------------------------" << std::endl;
     if (vsol->verbose > 2)
@@ -994,8 +1008,8 @@ void CoupledAdjoint::gradientwrtMesh(){
 Eigen::VectorXd CoupledAdjoint::runViscous()
 {
     int exitCode = vsol->run();
-    if (exitCode != 0)
-        std::cout << "vsol terminated with exit code " << exitCode << std::endl;
+    // if (exitCode != 0)
+    //     std::cout << "vsol terminated with exit code " << exitCode << std::endl;
     
     // Extract deltaStar
     std::vector<Eigen::VectorXd> dStar_p;