From 3b935a0e0aa431e54fe677eb42f84e71727639c6 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Sun, 29 Sep 2024 10:56:00 +0200
Subject: [PATCH] (feat) Added MDAO api for aerodynamics optimization
The api features 3 components: The solver, the geometry and the initial mesh. The api is tested in a aerodynamic shape case where the NACA0012 is constrained at 0deg AoA with cl=0.40 and the objective function to minimize is the drag wrt the max camber and the max camber position. The test case ends with a final shape close to the NACA4412 (camber = 35% & max camber pos = 41%)
---
blast/api/mdaAPI.py | 294 ++++++++++++++++++++++++++++++++++++
blast/mdao/naca0012_mdao.py | 201 ++++++++++++++++++++++++
2 files changed, 495 insertions(+)
create mode 100644 blast/api/mdaAPI.py
create mode 100644 blast/mdao/naca0012_mdao.py
diff --git a/blast/api/mdaAPI.py b/blast/api/mdaAPI.py
new file mode 100644
index 0000000..43fc618
--- /dev/null
+++ b/blast/api/mdaAPI.py
@@ -0,0 +1,294 @@
+#!/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.
+
+
+## BLASTER api for openMDAO
+# Paul Dechamps
+
+import numpy as np
+import openmdao.api as om
+import blast
+import blast.utils as viscUtils
+
+# Solver openMDAO object
+class BlasterSolver(om.ExplicitComponent):
+ """OpenMDAO component for the Blaster solver
+
+ Parameters
+ ----------
+ - cfg {dict} -- Configuration dictionary for the VII coupling
+ - icfg {dict} -- Configuration dictionary for the inviscid solver interface
+ - iSolverName {str} -- Inviscid solver name (default: 'DART')
+
+ Inputs
+ ------
+ - aoa {float} -- Angle of attack
+ - x_aero {np.array} -- Mesh coordinates
+
+ Outputs
+ -------
+ - cl {float} -- Lift coefficient
+ - cd {float} -- Drag coefficient
+ """
+
+ def __init__(self, cfg, icfg, iSolverName='DART'):
+ super().__init__()
+
+ self.nCpt = 0
+
+ if iSolverName == 'DART':
+ from blast.interfaces.dart.DartInterface import DartInterface as interface
+ else:
+ raise RuntimeError('Solver '+iSolverName+' currently not implemented')
+
+ for i in range(len(cfg['xtrF'])):
+ if cfg['xtrF'][i] is None:
+ cfg['xtrF'][i] = -1
+ self.vsol = blast.Driver(cfg['Re'], cfg['Minf'], cfg['CFL0'], cfg['nSections'], xtrF_top=cfg['xtrF'][0], xtrF_bot=cfg['xtrF'][1], _span=cfg['span'], _verbose=cfg['Verb'])
+
+ self.isol = interface(icfg, self.vsol, cfg)
+
+ # Coupler
+ import blast.coupler as blastCoupler
+ self.coupler = blastCoupler.Coupler(self.isol, self.vsol, _maxCouplIter=cfg['couplIter'], _couplTol=cfg['couplTol'], _iterPrint=cfg['iterPrint'], _resetInv=cfg['resetInv'])
+
+ # Adjoint
+ self.adjointSolver = blast.CoupledAdjoint(self.isol.adjointSolver, self.vsol)
+
+ # Single coupler run to initialize viscous sections (which depend on size of upper/lower sides)
+ self.coupler.run(write=False)
+ self.coupler.reset()
+ 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)
+
+ def setup_partials(self):
+ self.declare_partials('cl', 'aoa', method='exact')
+ self.declare_partials('cd', 'aoa', method='exact')
+ self.declare_partials('cl', 'x_aero', method='exact')
+ self.declare_partials('cd', 'x_aero', method='exact')
+
+
+ def compute_partials(self, inputs, partials, discrete_inputs=None):
+ """Compute partial derivatives using the adjoint solver
+ """
+ # Run adjoint solver
+ self.isol.iobj['adj'].run()
+ self.adjointSolver.reset()
+ self.adjointSolver.run()
+
+ # AoA derivatives
+ partials['cl', 'aoa'] = self.adjointSolver.tdCl_AoA
+ partials['cd', 'aoa'] = self.adjointSolver.tdCd_AoA
+
+ # 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):
+ """Update the mesh, aoa and run the solver
+ """
+ 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()
+
+ # Run forward problem
+ self.coupler.run(write=False)
+ outputs['cl'] = self.isol.getCl()
+ outputs['cd'] = self.isol.getCd() + self.vsol.Cdf
+
+ # Write results
+ self.isol.writeCp(sfx='_'+str(self.nCpt))
+ viscUtils.getSolution(self.isol.sec, write=True, toW=['x', 'y', 'cf'], sfx='_'+str(self.nCpt))
+ self.nCpt += 1
+
+ def getBoundaryMesh(self):
+ """Initial surface mesh coordinates
+ """
+ return BlasterMesh(nDim=self.isol.getnDim(), bnd=self.isol.iobj['bnd'])
+
+ def getBoundaryGeometry(self):
+ """Airfoil geometry
+ """
+ return BlasterGeometry(nDim=self.isol.getnDim(), nNodes=self.isol.iobj['bnd'].nodes.size(), nodeRows=self.isol.vBnd[0][0].connectListNodes)
+
+# Surface mesh
+class BlasterMesh(om.IndepVarComp):
+ """Initial 2D surface mesh coordinates
+
+ Parameters
+ ----------
+ - nDim {int} -- Number of dimensions
+ - bnd {object} -- Boundary object containing the mesh
+
+ Outputs
+ -------
+ - x_aero0 {np.array} -- Initial mesh coordinates
+ """
+
+ def initialize(self):
+ self.options.declare('nDim')
+ self.options.declare('bnd', recordable=False)
+
+ def setup(self):
+ bnd = self.options['bnd']
+ nDim = self.options['nDim']
+ x_aero0 = np.zeros(bnd.nodes.size() * nDim)
+ for inod, n in enumerate(bnd.nodes):
+ for idim in range(nDim):
+ x_aero0[inod*nDim + idim] = n.pos[idim]
+ self.add_output('x_aero0', val=x_aero0)
+
+class BlasterGeometry(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.
+ Parameters become tau (thickness), eps (max camber) and p (max camber position).
+
+ Parameters
+ ----------
+ - nDim {int} -- Number of dimensions
+ - nNodes {int} -- Number of nodes
+ - nodeRows {np.array} -- Nodal indexes
+
+ Inputs
+ ------
+ - x_aero0 {np.array} -- Initial mesh
+ - tau {float} -- Airfoil thickness
+ - eps {float} -- Airfoil max camber
+ - p {float} -- Airfoil max camber position
+
+ Outputs
+ -------
+ - x_aero_out {np.array} -- Airfoil coordinates
+ """
+
+ def initialize(self):
+ self.options.declare('nDim')
+ self.options.declare('nNodes')
+ self.options.declare('nodeRows')
+
+ 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('tau', val = 0.12, desc='Airfoil thickness')
+ self.add_input('eps', val = 0., desc='Airfoil max camber')
+ self.add_input('p', val = 0., desc='Airfoil max camber position')
+ self.add_output('x_aero_out', val = np.zeros(nNodes*nDim))
+
+ def setup_partials(self):
+ self.declare_partials('x_aero_out', 'tau', method='fd')
+ self.declare_partials('x_aero_out', 'eps', method='fd')
+ self.declare_partials('x_aero_out', 'p', method='fd')
+
+ def compute(self, inputs, outputs):
+ """Compute airfoil coordinates
+ Points are sorted in seilig format to identify upper and lower surfaces.
+ Then the NACA 4-digit airfoil is generated and the points are sorted back
+ to the original mesh format."""
+
+ x_aero0 = inputs['x_aero0']
+ tau = inputs['tau'][0]
+ eps = inputs['eps'][0]
+ p = inputs['p'][0]
+ nDim = self.options['nDim']
+
+ # Retreive selig coordinates
+ x_seilig = np.zeros(len(self.options['nodeRows']))
+ for i, val in enumerate(self.options['nodeRows']):
+ x_seilig[i] = x_aero0[val*nDim]
+ # Generate airfoil
+ # Format tau eps and p with 2 digits after coma
+ print(f"tau: {tau:.10f}, eps: {eps:.10f}, p: {p:.10f}")
+ xb, yb = self.__generateAirfoil(x_seilig, tau, eps, p)
+ # Go back to normal mesh
+ x_aero_out = np.zeros(len(x_aero0))
+ for i, val in enumerate(self.options['nodeRows']):
+ x_aero_out[val*nDim] = xb[i]
+ x_aero_out[val*nDim+1] = yb[i]
+ outputs['x_aero_out'] = x_aero_out
+
+ def __generateAirfoil(self, x, tau, eps, p):
+ """Generate NACA 4-digit airfoil from x coordinates points
+
+ Arguments
+ ---------
+ - x {np.array} -- x coordinates sorted in seilig format
+ - tau {float} -- Airfoil thickness
+ - eps {float} -- Airfoil max camber
+ - p {float} -- Airfoil max camber position
+
+ Returns
+ -------
+ - tuple -- x and y coordinates
+ """
+
+ T = 10 * tau * (
+ 0.2969 * np.sqrt(x)
+ - 0.1260 * (x)
+ - 0.3537 * (x)**2
+ + 0.2843 * (x)**3
+ - 0.1015 * (x)**4) # Thickness
+
+ YBar = np.zeros(len(T))
+ dYBar_dx = np.zeros(len(T))
+ Y = np.zeros(len(T))
+ xB = np.zeros(len(T))
+ yB = np.zeros(len(T))
+
+ for i, xi in enumerate(x):
+ # Camber
+ if xi < p:
+ YBar[i] = ((eps * xi)/(p**2)) * (2 * p - xi)
+ dYBar_dx[i] = (2 * eps)/(p**2) * (p - xi)
+ elif xi >= p:
+ YBar[i] = (eps * (1 - xi)/((1-p)**2)) * (1 + xi - 2 * p)
+ dYBar_dx[i] = (2 * eps)/((1 - p)**2) * (p - xi)
+
+ # Final x & y coordinates
+ if i >= len(x)/2 - 1:
+ # Lower surface
+ 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]))
+ else:
+ # Upper surface
+ 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)
diff --git a/blast/mdao/naca0012_mdao.py b/blast/mdao/naca0012_mdao.py
new file mode 100644
index 0000000..ec34a6e
--- /dev/null
+++ b/blast/mdao/naca0012_mdao.py
@@ -0,0 +1,201 @@
+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' : 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.2, # freestream Mach number
+ 'AoA' : 2., # 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.2, # 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
+ }
+
+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.getBoundaryGeometry())
+ 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=0.0)
+ self.connect('aoa', 'solver.aoa')
+ #self.add_design_var('aoa', lower=0.0, upper=3.0)
+
+ self.dvs.add_output('tau', val=0.12)
+ self.connect('tau', 'geometry.tau')
+ #self.add_design_var('tau', lower=0.10, upper=0.2)
+
+ 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.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_constraint('solver.cl', equals=0.4)
+
+ #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()
+
+ # 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()
--
GitLab