From 046b10187b6cbae745c53c532b664fcabd5ae7d5 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Wed, 13 Nov 2024 17:40:06 +0100
Subject: [PATCH] (tests) Updated aero optimization cases

Refactor api + updated tests with mode updater and added scalers
---
 blast/api/mdaAPI.py          |  4 ----
 blast/mdao/opti_NACA.py      | 17 ++++++++++-------
 blast/mdao/opti_subsonic.py  | 37 ++++++++++++++++++++++++------------
 blast/mdao/opti_transonic.py | 16 ++++++++--------
 4 files changed, 43 insertions(+), 31 deletions(-)

diff --git a/blast/api/mdaAPI.py b/blast/api/mdaAPI.py
index 8a4cae9..595572e 100644
--- a/blast/api/mdaAPI.py
+++ b/blast/api/mdaAPI.py
@@ -241,16 +241,12 @@ class BlasterModeUpdater(om.ExplicitComponent):
 
         # 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']
diff --git a/blast/mdao/opti_NACA.py b/blast/mdao/opti_NACA.py
index 4c9daeb..e14e0f2 100644
--- a/blast/mdao/opti_NACA.py
+++ b/blast/mdao/opti_NACA.py
@@ -7,6 +7,7 @@ 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
 
@@ -36,8 +37,8 @@ def cfgInviscid(nthrds, verb):
     'dbc' : False,
     'Upstream' : 'upstream',
     # Freestream
-    'M_inf' : 0.2, # freestream Mach number
-    'AoA' : 2., # freestream angle of attack
+    '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
@@ -84,6 +85,7 @@ class groupMDA(om.Group):
 
         # Design variables box
         self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*'])
+        self.add_subsystem('mod', BlasterModeUpdater(tolerance=10., force_mode='visc')) # Mode updater
         self.add_subsystem('mesh', blaster.getBoundaryMesh())
         self.add_subsystem('geometry', blaster.getBoundaryNACAGeometry())
         self.add_subsystem('solver', blaster)
@@ -94,6 +96,7 @@ class groupMDA(om.Group):
         # 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')
 
         #nonlinear_solver = om.NewtonSolver(solve_subsystems=False)
         self.linear_solver = om.DirectSolver()
@@ -102,7 +105,7 @@ class groupMDA(om.Group):
         # 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.add_design_var('aoa', lower=-3.0, upper=3.0, scaler=1)
 
         self.dvs.add_output('tau', val=0.12)
         self.connect('tau', 'geometry.tau')
@@ -110,17 +113,17 @@ class groupMDA(om.Group):
 
         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.add_design_var('eps', lower=0.0, upper=0.05, scaler=10)
 
         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_design_var('p', lower=0.35, upper=0.45, scaler=10)
 
-        self.add_constraint('solver.cl', equals=0.4)
+        self.add_constraint('solver.cl', equals=0.4, scaler=10)
 
         #self.connect('solver.cl', 'cl')
         self.connect('solver.cd', 'cd')
-        self.add_objective('obj')
+        self.add_objective('obj', scaler=1000)
 
 def main():
     # Timer.
diff --git a/blast/mdao/opti_subsonic.py b/blast/mdao/opti_subsonic.py
index bbe4496..2228498 100644
--- a/blast/mdao/opti_subsonic.py
+++ b/blast/mdao/opti_subsonic.py
@@ -7,6 +7,7 @@ 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
 
@@ -23,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' : 25, 'msTe' : 0.0075, 'msLe' : 0.001}, # parameters for input file model
+    'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF' : 25, 'msTe' : 0.01, 'msLe' : 0.005}, # parameters for input file model
     'Dim' : 2, # problem dimension
     'Format' : 'gmsh', # save format (vtk or gmsh)
     # Markers
@@ -37,7 +38,7 @@ def cfgInviscid(nthrds, verb):
     'Upstream' : 'upstream',
     # Freestream
     'M_inf' : 0.1, # freestream Mach number
-    'AoA' : 1., # freestream angle of attack
+    'AoA' : 0., # freestream angle of attack
     # Geometry
     'S_ref' : 1., # reference surface length
     'c_ref' : 1., # reference chord length
@@ -84,10 +85,17 @@ class groupMDA(om.Group):
 
         # Design variables box
         self.add_subsystem('dvs', om.IndepVarComp(), promotes=['*'])
+        self.add_subsystem('mod', BlasterModeUpdater(tolerance=10., force_mode='visc')) # Mode updater
         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]',
+        self.add_subsystem('constraint_comp_0', om.ExecComp('constraint_var_0 = lower_shape[0] + upper_shape[0]',
+                                                          lower_shape=np.zeros(4), upper_shape=np.zeros(4)), promotes=['*'])
+        self.add_subsystem('constraint_comp_1', om.ExecComp('constraint_var_1 = lower_shape[1] + upper_shape[1]',
+                                                          lower_shape=np.zeros(4), upper_shape=np.zeros(4)), promotes=['*'])
+        self.add_subsystem('constraint_comp_2', om.ExecComp('constraint_var_2 = lower_shape[2] + upper_shape[2]',
+                                                          lower_shape=np.zeros(4), upper_shape=np.zeros(4)), promotes=['*'])
+        self.add_subsystem('constraint_comp_3', om.ExecComp('constraint_var_3 = lower_shape[3] + upper_shape[3]',
                                                           lower_shape=np.zeros(4), upper_shape=np.zeros(4)), promotes=['*'])
         self.add_subsystem('solver', blaster)
 
@@ -97,23 +105,25 @@ class groupMDA(om.Group):
         # 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=1.0)
+        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('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.add_design_var('upper_shape', lower=-0.1, upper=0.4, scaler=10)
 
         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.add_design_var('lower_shape', lower=-0.4, upper=0.1, scaler=10)
 
         self.dvs.add_output('class_shape_n1', val=0.5)
         self.connect('class_shape_n1', 'geometry.class_shape_n1')
@@ -125,14 +135,17 @@ class groupMDA(om.Group):
 
         # 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.add_constraint('solver.cl', equals=0., scaler=10)
+        #self.add_constraint('geometry.thickness', upper=.1)
+        self.add_constraint('geometry.LEradius', equals=0.0224, scaler=100)
+        self.add_constraint('constraint_var_0', equals=0., linear=True)
+        self.add_constraint('constraint_var_1', equals=0., linear=True)
+        self.add_constraint('constraint_var_2', equals=0., linear=True)
+        self.add_constraint('constraint_var_3', 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.
@@ -144,7 +157,7 @@ def main():
 
     prob.driver = om.ScipyOptimizeDriver()
     prob.driver.options['optimizer'] = 'SLSQP'
-    prob.driver.options['tol'] = 1e-3
+    prob.driver.options['tol'] = 1e-5
 
     recorder = om.SqliteRecorder('reports/BlasterSolver_file.sql')
 
diff --git a/blast/mdao/opti_transonic.py b/blast/mdao/opti_transonic.py
index 5cbdeca..0d50a1a 100644
--- a/blast/mdao/opti_transonic.py
+++ b/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.002}, # parameters for input file model
+    'Pars' : {'xLgt' : 50, 'yLgt' : 50, 'msF' : 20, 'msTe' : 0.01, 'msLe' : 0.007}, # parameters for input file model
     'Dim' : 2, # problem dimension
     'Format' : 'gmsh', # save format (vtk or gmsh)
     # Markers
@@ -37,8 +37,8 @@ def cfgInviscid(nthrds, verb):
     'dbc' : False,
     'Upstream' : 'upstream',
     # Freestream
-    'M_inf' : 0.1, # freestream Mach number
-    'AoA' : 0., # freestream angle of attack
+    'M_inf' : 0.78, # freestream Mach number
+    'AoA' : 2., # freestream angle of attack
     # Geometry
     'S_ref' : 1., # reference surface length
     'c_ref' : 1., # reference chord length
@@ -52,17 +52,17 @@ def cfgInviscid(nthrds, verb):
     '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
+    'Max_it' : 50 # 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)
+        '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-5,          # Tolerance of the VII methodology
+        '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
@@ -77,7 +77,7 @@ class groupMDA(om.Group):
         args = parseargs()
         icfg = cfgInviscid(args.k, args.verb)
         vcfg = cfgBlast(args.verb)
-        gr = 1.055
+        gr = 1.1
         if gr == 1.0:
             icfg['Pars']['msF'] = icfg['Pars']['msLe']
         else:
@@ -134,7 +134,7 @@ class groupMDA(om.Group):
 
         # Constraints
         # Cl constraint
-        self.add_constraint('solver.cl', equals=0.7)
+        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('constraint_var', equals=0., linear=True)
-- 
GitLab