From 31238be0c84fb7c983b3b3634c6f8e2a0fa137f8 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Mon, 27 Jan 2025 19:19:37 +0100
Subject: [PATCH] (feat) Refactor optimzation problem setup

---
 blast/api/mda_api.py          | 109 ++++++++++++++++++-
 blast/mdao/opti_NACA2.py      | 196 ++++++++++++++++++++++++++++++++++
 blast/mdao/opti_transonic.py  |   2 +-
 blast/mdao/opti_transonic2.py | 179 +++++++++++++++++++++++++++++++
 4 files changed, 483 insertions(+), 3 deletions(-)
 create mode 100644 blast/mdao/opti_NACA2.py
 create mode 100644 blast/mdao/opti_transonic2.py

diff --git a/blast/api/mda_api.py b/blast/api/mda_api.py
index 55b5bab..dbff35b 100644
--- a/blast/api/mda_api.py
+++ b/blast/api/mda_api.py
@@ -24,7 +24,6 @@ from fwk.coloring import ccolors
 import openmdao.api as om
 import blast
 import blast.blUtils as vutils
-from pygeo import DVGeometryCST
 from scipy.optimize import least_squares
 import fwk
 
@@ -79,7 +78,7 @@ class BlasterSolver(om.ExplicitComponent):
         self.coupler.run(write=False)
         self.__write(sfx='baseline')
         self.coupler.reset()
-        print('object created')
+        print('BLASTER objects created')
     
     def setup(self):
         self.add_input('aoa', val=0.0, desc='Flow angle of attack')
@@ -205,6 +204,7 @@ class BlasterSolver(om.ExplicitComponent):
     def getBoundaryFFDGeometry(self):
         """Airfoil FFD geometry
         """
+        from pygeo import DVGeometryCST
         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")
@@ -498,3 +498,108 @@ class BlasterFFDGeometry(om.ExplicitComponent):
         Ri = calc_R(*center)
         R = Ri.mean()
         return center, R
+
+class groupMDA(om.Group):
+    def __init__(self, icfg, vcfg, ocfg):
+        super().__init__()
+        self.icfg = icfg
+        self.vcfg = vcfg
+        self.ocfg = ocfg
+
+        if 'aoa' not in self.ocfg['initialValues']:
+            raise ValueError('Initial value for aoa is missing')
+        if self.ocfg['geometryType'] == 'NACA':
+            if 'tau' not in self.ocfg['initialValues']:
+                raise ValueError('Initial value for tau is missing')
+            if 'eps' not in self.ocfg['initialValues']:
+                raise ValueError('Initial value for eps is missing')
+            if 'p' not in self.ocfg['initialValues']:
+                raise ValueError('Initial value for p is missing')
+            for key in self.ocfg['initialValues']:
+                if key not in ['aoa', 'tau', 'eps', 'p']:
+                    raise KeyError('Invalid key in initialValues', key)
+        elif self.ocfg['geometryType'] == 'FFD':
+            if 'upper_shape' not in self.ocfg['initialValues']:
+                raise ValueError('Initial value for upper_shape is missing')
+            if 'lower_shape' not in self.ocfg['initialValues']:
+                raise ValueError('Initial value for lower_shape is missing')
+            if 'class_shape_n1' not in self.ocfg['initialValues']:
+                raise ValueError('Initial value for class_shape_n1 is missing')
+            if 'class_shape_n2' not in self.ocfg['initialValues']:
+                raise ValueError('Initial value for class_shape_n2 is missing')
+            for key in self.ocfg['initialValues']:
+                if key not in ['aoa', 'upper_shape', 'lower_shape', 'class_shape_n1', 'class_shape_n2']:
+                    raise KeyError('Invalid key in initialValues', key)
+
+    def setup(self):
+        # Add the MDAO components
+        blaster = BlasterSolver(self.vcfg, self.icfg)
+
+        print('Optimization setup')
+        print('------------------')
+        # Design variables box
+        self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*'])
+        self.add_subsystem('mod', BlasterModeUpdater(tolerance=10., force_mode=self.ocfg['mode'])) # Mode updater
+        self.add_subsystem('mesh', blaster.getBoundaryMesh())
+        # Geometry component
+        if self.ocfg['geometryType'] == 'NACA':
+            print('Adding NACA geometry')
+            self.add_subsystem('geometry', blaster.getBoundaryNACAGeometry())
+        elif self.ocfg['geometryType'] == 'FFD':
+            print('Adding FFD geometry')
+            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(len(self.ocfg['initialValues']['upper_shape']['value'])), upper_shape=np.zeros(len(self.ocfg['initialValues']['lower_shape']['value']))), promotes=['*'])
+        else:
+            raise ValueError('Invalid geometryType', self.ocfg['geometryType'])
+        # Solver
+        self.add_subsystem('solver', blaster)
+        # Objectives
+        for obj in self.ocfg['objectives']:
+            print('Adding objective', obj['name'])
+            self.add_subsystem('obj_cmp_'+obj['name'], obj['object'], promotes=obj['promotions'])
+
+        # Connect meshes
+        self.connect('mesh.x_aero0', 'geometry.x_aero0')
+        self.connect('geometry.x_aero_out', 'solver.x_aero')
+        self.connect('mod.visc_mode', 'solver.visc_mode')
+
+        if 'nonlinear_solver' in self.ocfg:
+            self.nonlinear_solver = self.ocfg['nonlinearSolver']
+        if 'linear_solver' in self.ocfg:
+            self.linear_solver = self.ocfg['linearSolver']
+        print('')
+    
+    def configure(self):
+        print('Optimization config')
+        print('------------------')
+        for variable, setup in self.ocfg['initialValues'].items():
+            print(f'- dvs output {variable} linked to {setup["link"]} with value {setup["value"]}')
+            self.dvs.add_output(variable, val=setup['value'])
+            self.connect(variable, setup['link']+'.'+variable)
+            if variable in self.ocfg['designVariables']:
+                design_config = self.ocfg['designVariables'][variable]
+                print(f'- design variable {variable} with bounds {design_config["lower_bound"]} {design_config["upper_bound"]} and scaler {design_config["scaler"]}')
+                self.add_design_var(variable, lower=design_config['lower_bound'], upper=design_config['upper_bound'], scaler=design_config['scaler'])
+        
+        for const in self.ocfg['constraints']:
+            if const['type'] == 'equality':
+                print('- equality constraint:', const['name'], '=', const['value'])
+                self.add_constraint(const['name'], equals=const['value'], scaler=const['scaler'])
+            elif const['type'] == 'inequality':
+                print('- inequality constraint:', const['name'], end=' ')
+                if 'lower' in const:
+                    print('-- lower bound:', const['lower'])
+                    self.add_constraint(const['name'], lower=const['lower'], scaler=const['scaler'])
+                elif 'upper' in const:
+                    print('-- upper bound:', const['upper'])
+                    self.add_constraint(const['name'], upper=const['upper'], scaler=const['scaler'])
+        if self.ocfg['geometryType'] == 'FFD':
+            self.add_constraint('constraint_var', equals=0.0, linear=True)
+
+        for obj in self.ocfg['objectives']:
+            print('- objective:', obj['name'], 'linked to', obj['link'])
+            self.connect(obj['link']+'.'+obj['name'], obj['name'])
+            self.add_objective('obj', scaler=obj['scaler'])
+        print('')
diff --git a/blast/mdao/opti_NACA2.py b/blast/mdao/opti_NACA2.py
new file mode 100644
index 0000000..8666f19
--- /dev/null
+++ b/blast/mdao/opti_NACA2.py
@@ -0,0 +1,196 @@
+#!/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.
+
+
+# @author Paul Dechamps
+# @date 2024
+# Optimization case using NACA parametrization
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+from blast.api.mda_api import groupMDA
+
+import numpy as np
+import openmdao.api as om
+
+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.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' : '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.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-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
+    }
+
+def cfgOpti():
+    return{
+        'geometryType': 'NACA',
+        'mode': 'visc',
+        'linearSolver': om.DirectSolver(),
+        'objectives': [{'name': 'cd',
+                        'link': 'solver',
+                        'object': om.ExecComp('obj = cd'),
+                        'promotions': ['obj', 'cd'],
+                        'scaler': 1000}
+                        ],
+        'initialValues': {'aoa':{'link':'solver', 'value': 0.0},
+                            'tau':{'link':'geometry', 'value': 0.12},
+                            'eps':{'link':'geometry', 'value': 0.0},
+                            'p':{'link': 'geometry', 'value': 0.4}
+        },
+        'designVariables': {'aoa':{'lower_bound': -3.0, 'upper_bound': 3.0, 'scaler': 1},
+                            'eps':{'lower_bound': 0.0, 'upper_bound': 0.05, 'scaler': 10},
+                            'p':{'lower_bound': 0.35, 'upper_bound': 0.45, 'scaler': 10}
+        },
+        'constraints': [{'name': 'solver.cl', 'type': 'equality', 'value': 0.4, 'scaler': 10}]
+    }
+
+def main():
+    # Timer.
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    args = parseargs()
+    icfg = cfgInviscid(args.k, args.verb)
+    vcfg = cfgBlast(args.verb)
+    ocfg = cfgOpti()
+
+    prob = om.Problem()
+    prob.model = groupMDA(icfg, vcfg, ocfg)
+
+    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()
diff --git a/blast/mdao/opti_transonic.py b/blast/mdao/opti_transonic.py
index 7ca7ffe..7b74ed3 100644
--- a/blast/mdao/opti_transonic.py
+++ b/blast/mdao/opti_transonic.py
@@ -63,7 +63,7 @@ def cfgInviscid(nthrds, verb):
     '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)
+    '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
diff --git a/blast/mdao/opti_transonic2.py b/blast/mdao/opti_transonic2.py
new file mode 100644
index 0000000..68c695f
--- /dev/null
+++ b/blast/mdao/opti_transonic2.py
@@ -0,0 +1,179 @@
+#!/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.
+
+
+# @author Paul Dechamps
+# @date 2024
+# Optimization for a transonic case
+
+from fwk.wutils import parseargs
+import fwk
+from fwk.testing import *
+from fwk.coloring import ccolors
+
+from blast.api.mda_api import groupMDA
+
+import numpy as np
+import openmdao.api as om
+
+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.008}, # 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
+    'Upstream' : 'upstream',
+    # Freestream
+    'M_inf' : 0.78, # 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' : 50 # solver maximum number of iterations
+    }
+
+def cfgBlast(verb):
+    return {
+        'Re' : 1e8,                 # Freestream Reynolds number
+        'Minf' : 0.78,               # 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-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' : [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
+    }
+
+def cfgOpti():
+    return{
+        'geometryType': 'FFD',
+        'mode': 'visc',
+        'linearSolver': om.DirectSolver(),
+        'objectives': [{'name': 'cd',
+                        'link': 'solver',
+                        'object': om.ExecComp('obj = cd'),
+                        'promotions': ['obj', 'cd'],
+                        'scaler': 1000}
+                        ],
+        'initialValues': {'aoa':{'link':'solver', 'value': 0.0},
+                            'upper_shape':{'link':'geometry', 'value': [0.17553542, 0.14723497, 0.14531762, 0.13826955]},
+                            'lower_shape':{'link':'geometry', 'value': [-0.16935506, -0.15405552, -0.13912672, -0.14168947]},
+                            'class_shape_n1':{'link': 'geometry', 'value': 0.5},
+                            'class_shape_n2':{'link': 'geometry', 'value': 1.0}
+        },
+        'designVariables':  {'aoa':{'lower_bound': -5.0, 'upper_bound': 5.0, 'scaler': 1.0},
+                            'upper_shape':{'lower_bound': -0.15, 'upper_bound': 0.25, 'scaler': 10},
+                            'lower_shape':{'lower_bound': -0.2, 'upper_bound': 0.2, 'scaler': 10}
+        },
+        'constraints': [{'name': 'solver.cl', 'type': 'equality', 'value': 0.4, 'scaler': 10, 'linear': False},
+                        {'name': 'geometry.thickness', 'type': 'inequality', 'lower': 0.12, 'scaler':10, 'upper':None, 'linear': True},
+                        {'name': 'geometry.LEradius', 'type': 'equality', 'value': 0.0224, 'scaler':100, 'linear': True}]
+    }
+
+def main():
+    # Timer.
+    tms = fwk.Timers()
+    tms['total'].start()
+
+    args = parseargs()
+    icfg = cfgInviscid(args.k, args.verb)
+    vcfg = cfgBlast(args.verb)
+    ocfg = cfgOpti()
+
+    prob = om.Problem()
+    prob.model = groupMDA(icfg, vcfg, ocfg)
+
+    prob.driver = om.ScipyOptimizeDriver()
+    prob.driver.options['optimizer'] = 'SLSQP'
+    prob.driver.options['tol'] = 1e-4
+
+    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()
+
+    # problem timers
+    print(prob.model.solver.tms)
+
+    # Case reader
+    cr = om.CaseReader(prob.get_reports_dir() / "../BlasterSolver_file.sql")
+    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('aoa:', prob['aoa'])
+    print('cl:', prob['solver.cl'])
+    print('cd:', prob['solver.cd'])
+    print('obj:', prob['obj'])
+
+    # eof
+    tms['total'].stop()
+    print(tms)
+
+if __name__ == "__main__":
+    main()
-- 
GitLab