(fix) Added possibility to save all steps in optimization + refactor transonic optimization

blast/api/mdaAPI.py

@@ -49,9 +49,10 @@ class BlasterSolver(om.ExplicitComponent):
     - cd {float} -- Drag coefficient
     """
 
-    def __init__(self, cfg, icfg, iSolverName='DART'):
+    def __init__(self, cfg, icfg, iSolverName='DART', saveAll=False):
         super().__init__()
 
+        self._saveAll = saveAll
         self.nCpt = 0
         self.tms = fwk.Timers()
 
@@ -76,6 +77,7 @@ class BlasterSolver(om.ExplicitComponent):
 
         # Single coupler run to initialize viscous sections (which depend on size of upper/lower sides)
         self.coupler.run(write=False)
+        self.__write(sfx='baseline')
         self.coupler.reset()
         print('object created')
     
@@ -172,8 +174,8 @@ class BlasterSolver(om.ExplicitComponent):
             raise ValueError('omBlasterSolver::compute - Invalid mode', inputs['visc_mode'][0])
         self.tms['fwd'].stop()
 
-        # write solution
-        self.__write()
+        # write solutions
+        self.__write(str(self.nCpt) if self._saveAll else 'opt')
 
         outputs['cl'] = self.isol.getCl()
         outputs['cd'] = self.isol.getCd()
@@ -181,14 +183,8 @@ class BlasterSolver(om.ExplicitComponent):
             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'
+    def __write(self, sfx='0'):
+        sfx = '_'+sfx
         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

blast/mdao/opti_transonic.py

@@ -24,7 +24,7 @@ def cfgInviscid(nthrds, verb):
     '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.007}, # parameters for input file 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
@@ -38,7 +38,7 @@ def cfgInviscid(nthrds, verb):
     'Upstream' : 'upstream',
     # Freestream
     'M_inf' : 0.78, # freestream Mach number
-    'AoA' : 2., # freestream angle of attack
+    'AoA' : 0., # freestream angle of attack
     # Geometry
     'S_ref' : 1., # reference surface length
     'c_ref' : 1., # reference chord length
@@ -62,7 +62,7 @@ def cfgBlast(verb):
         'CFL0' : 1.,                # Inital CFL number of the calculation
         'Verb': verb,               # Verbosity level of the solver
         'couplIter': 75,            # Maximum number of iterations
-        'couplTol' : 1e-3,          # Tolerance of the VII methodology
+        '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
@@ -85,7 +85,7 @@ class groupMDA(om.Group):
             icfg['Pars']['msF'] = icfg['Pars']['msLe']*gr**(n-1)
 
         # Add the MDAO components
-        blaster = BlasterSolver(vcfg, icfg)
+        blaster = BlasterSolver(vcfg, icfg, iSolverName='DART', saveAll=False)
 
         self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*']) # Design variables box
         self.add_subsystem('mod', BlasterModeUpdater(tolerance=10., force_mode='visc')) # Mode updater
@@ -110,17 +110,17 @@ class groupMDA(om.Group):
         # 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.add_design_var('aoa', lower=-5, upper=5)
 
         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.add_design_var('upper_shape', lower=-0.15, upper=0.25, scaler=10)
 
         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.add_design_var('lower_shape', lower=-0.2, upper=0.2, scaler=10)
 
         self.dvs.add_output('class_shape_n1', val=0.5)
         self.connect('class_shape_n1', 'geometry.class_shape_n1')
@@ -130,18 +130,14 @@ class groupMDA(om.Group):
         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.4)
-        self.add_constraint('geometry.thickness', lower=0.12)
-        self.add_constraint('geometry.LEradius', equals=0.0224)
+        self.add_constraint('solver.cl', equals=0.4, scaler=10)
+        self.add_constraint('geometry.thickness', lower=0.12, linear=True, scaler=10)
+        self.add_constraint('geometry.LEradius', equals=0.0224, linear=True, scaler=100)
         self.add_constraint('constraint_var', equals=0., linear=True)
 
-        #self.connect('solver.cl', 'cl')
         self.connect('solver.cd', 'cd')
-        self.add_objective('obj')
+        self.add_objective('obj', scaler=1000)
 
 def main():
     # Timer.
@@ -153,7 +149,7 @@ def main():
 
     prob.driver = om.ScipyOptimizeDriver()
     prob.driver.options['optimizer'] = 'SLSQP'
-    prob.driver.options['tol'] = 1e-3
+    prob.driver.options['tol'] = 1e-4
 
     recorder = om.SqliteRecorder('reports/BlasterSolver_file.sql')
 
@@ -163,16 +159,16 @@ def main():
     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')
+    # print('check')
+    # prob.check_partials(includes=['cl', 'cd'], compact_print=True)
+    # print('done')
     prob.run_driver()
 
+    # problem timers
     print(prob.model.solver.tms)
 
-    # Instantiate your CaseReader
+    # Case reader
     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
@@ -190,37 +186,14 @@ def main():
         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()
-
+    # eof
+    tms['total'].stop()
+    print(tms)
 
 if __name__ == "__main__":
     main()