diff --git a/blast/api/mdaAPI.py b/blast/api/mdaAPI.py
index 43fc618cf91f0024418d3b763a26f61f3c9d11db..8a4cae98f25e36d9aea21f5cf2499de41a5abe1e 100644
--- a/blast/api/mdaAPI.py
+++ b/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
diff --git a/blast/mdao/naca0012_mdao.py b/blast/mdao/opti_NACA.py
similarity index 99%
rename from blast/mdao/naca0012_mdao.py
rename to blast/mdao/opti_NACA.py
index ec34a6ead36e2ad5799ead7bb2584b3b4f2385ba..4c9daebdf11d35b78ff94d6156a86fd4681bae72 100644
--- a/blast/mdao/naca0012_mdao.py
+++ b/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
diff --git a/blast/mdao/opti_subsonic.py b/blast/mdao/opti_subsonic.py
new file mode 100644
index 0000000000000000000000000000000000000000..bbe44968afb258523e5be46855d5303f40172497
--- /dev/null
+++ b/blast/mdao/opti_subsonic.py
@@ -0,0 +1,215 @@
+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()
diff --git a/blast/mdao/opti_transonic.py b/blast/mdao/opti_transonic.py
new file mode 100644
index 0000000000000000000000000000000000000000..5cbdeca570d336aa4ecd0760a0f99bf4b2951b1f
--- /dev/null
+++ b/blast/mdao/opti_transonic.py
@@ -0,0 +1,226 @@
+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()
diff --git a/blast/src/DCoupledAdjoint.cpp b/blast/src/DCoupledAdjoint.cpp
index 4cb1449ab6b6acacdaaa12f225c187019582ed4d..90d440a9dc75ed595ab5075f764d4cf62b1f78c4 100644
--- a/blast/src/DCoupledAdjoint.cpp
+++ b/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;
diff --git a/logo/logo.png b/logo/logo.png
index 52167245921b0e54b21c69be15c15829eac13dd7..ab93625029c9a6c45c6c48169588bbd72e9a538d 100644
Binary files a/logo/logo.png and b/logo/logo.png differ