Aerodynamic Optimization
Introduction
We will now demonstrate how to optimize the aerodynamic shape of a wing. We will be combining aspects of all of the following sections: Analysis with ADflow, pyOptSparse, and Geometric Parametrization. Again reference a high level explanation of how these fit together - aero problem, optimization problem, geometry. The optimization problem is defined as
Files
Navigate to the directory opt/aero in your tutorial folder.
Copy the following files to this directory:
cp ../ffd/ffd.xyz .
cp ../../aero/meshing/volume/wing_vol.cgns .
Create the following empty runscript in the current directory:
aero_opt.py
Dissecting the aerodynamic optimization script
Open the file aero_opt.py in your favorite text editor.
Then copy the code from each of the following sections into this file.
Import libraries
import os
import argparse
import ast
import numpy as np
from mpi4py import MPI
from baseclasses import AeroProblem
from adflow import ADFLOW
from pygeo import DVGeometry, DVConstraints
from pyoptsparse import Optimization, OPT
from idwarp import USMesh
from multipoint import multiPointSparse
The multipoint library is the only new library to include in this script.
Adding command line arguments
This is a convenience feature that allows the user to pass in command line arguments to the script. Two options are provided:
specifying the output directory
specifying the optimizer to be used
# Use Python's built-in Argument parser to get commandline options
parser = argparse.ArgumentParser()
parser.add_argument("--output", type=str, default="output")
parser.add_argument("--opt", type=str, default="SLSQP", choices=["IPOPT", "SLSQP", "SNOPT"])
parser.add_argument("--gridFile", type=str, default="wing_vol.cgns")
parser.add_argument("--optOptions", type=ast.literal_eval, default={}, help="additional optimizer options to be added")
args = parser.parse_args()
Creating processor sets
The multiPointSparse class allows us to allocate sets of processors for different analyses. This can be helpful if we want to consider multiple design points, each with a different set of flow conditions. In this case, we create a processor set for cruise cases, but we only add one point.
MP = multiPointSparse(MPI.COMM_WORLD)
MP.addProcessorSet("cruise", nMembers=1, memberSizes=MPI.COMM_WORLD.size)
comm, setComm, setFlags, groupFlags, ptID = MP.createCommunicators()
if not os.path.exists(args.output):
if comm.rank == 0:
os.mkdir(args.output)
If we want to add more points, we can increase the quantity nMembers.
We can choose the number of processors per point with the argument memberSizes.
We can add another processor set with another call to addProcessorSet.
The call createCommunicators returns information about the current processor’s set.
ADflow set-up
The set-up for adflow should look the same as for the aerodynamic analysis script. We add a single lift distribution with 200 sampling points.
aeroOptions = {
# I/O Parameters
"gridFile": args.gridFile,
"outputDirectory": args.output,
"monitorvariables": ["resrho", "cl", "cd"],
"writeTecplotSurfaceSolution": True,
# Physics Parameters
"equationType": "RANS",
# Solver Parameters
"smoother": "DADI",
"MGCycle": "sg",
"infchangecorrection": True,
# ANK Solver Parameters
"useANKSolver": True,
# NK Solver Parameters
"useNKSolver": True,
"nkswitchtol": 1e-6,
# Termination Criteria
"L2Convergence": 1e-10,
"L2ConvergenceCoarse": 1e-2,
"nCycles": 10000,
# Adjoint Parameters
"adjointL2Convergence": 1e-10,
}
# Create solver
CFDSolver = ADFLOW(options=aeroOptions, comm=comm)
CFDSolver.addLiftDistribution(150, "z")
CFDSolver.addSlices("z", np.linspace(0.1, 14, 10))
Set the AeroProblem
ap = AeroProblem(name="wing", alpha=1.5, mach=0.8, altitude=10000, areaRef=45.5, chordRef=3.25, evalFuncs=["cl", "cd"])
# Add angle of attack variable
ap.addDV("alpha", value=1.5, lower=0, upper=10.0, scale=0.1)
The only difference in setting up the AeroProblem is that now we add angle-of-attack as a design variable. Any of the quantities included in the initialization of the AeroProblem can be added as design variables.
Geometric parametrization
The set-up for DVGeometry should look very familiar (if not, see Geometric Parametrization).
We include twist and local variables in the optimization.
After setting up the DVGeometry instance we have to provide it to ADflow with the call setDVGeo.
# Create DVGeometry object
FFDFile = "ffd.xyz"
DVGeo = DVGeometry(FFDFile)
# Create reference axis
nRefAxPts = DVGeo.addRefAxis("wing", xFraction=0.25, alignIndex="k")
nTwist = nRefAxPts - 1
# Set up global design variables
def twist(val, geo):
for i in range(1, nRefAxPts):
geo.rot_z["wing"].coef[i] = val[i - 1]
DVGeo.addGlobalDV(dvName="twist", value=[0] * nTwist, func=twist, lower=-10, upper=10, scale=0.01)
# Set up local design variables
DVGeo.addLocalDV("local", lower=-0.5, upper=0.5, axis="y", scale=1)
# Add DVGeo object to CFD solver
CFDSolver.setDVGeo(DVGeo)
Geometric constraints
We can set up constraints on the geometry with the DVConstraints class, also found in the pyGeo module. There are several built-in constraint functions within the DVConstraints class, including thickness, surface area, volume, location, and general linear constraints. The majority of the constraints are defined based on a triangulated-surface representation of the wing obtained from ADflow.
Note
The triangulated surface is created by ADflow (or DAfoam) using the wall surfaces defined in the CFD volume mesh. The resolution is similar to the CFD surface mesh, and users do not need to provide this triangulated mesh themselves.
Optionally, this can also be defined with an external file, see the docstrings for setSurface().
This is useful if users want to have a different resolution on the triangulated surface (finer or coarser) compared to the CFD mesh, or if DVConstraints is being used without ADflow (or DAfoam).
The volume and thickness constraints are set up by creating a uniformly spaced 2D grid of points, which is then projected onto the upper and lower surface of a triangulated-surface representation of the wing.
The grid is defined by providing four corner points (using leList and teList) and by specifying the number of spanwise and chordwise points (using nSpan and nChord).
Note
These grid points are projected onto the triangulated surface along the normals of the ruled surface formed by these grid points.
Typically, leList and teList are given such that the two curves lie in a plane.
This ensures that the projection vectors are always exactly normal to this plane.
If the surface formed by leList and teList is not planar, issues can arise near the end of an open surface (i.e., the root of a wing) which can result in failing intersections.
By default, scaled=True for addVolumeConstraint() and addThicknessConstraints2D(), which means that the volume and thicknesses calculated will be scaled by the initial values (i.e., they will be normalized).
Therefore, lower=1.0 in this example means that the lower limits for these constraints are the initial values (i.e., if lower=0.5 then the lower limits would be half the initial volume and thicknesses).
Warning
The leList and teList points must lie completely inside the wing.
DVCon = DVConstraints()
DVCon.setDVGeo(DVGeo)
# Only ADflow has the getTriangulatedSurface Function
DVCon.setSurface(CFDSolver.getTriangulatedMeshSurface())
# Volume constraints
leList = [[0.01, 0, 0.001], [7.51, 0, 13.99]]
teList = [[4.99, 0, 0.001], [8.99, 0, 13.99]]
DVCon.addVolumeConstraint(leList, teList, nSpan=20, nChord=20, lower=1.0, scaled=True)
# Thickness constraints
DVCon.addThicknessConstraints2D(leList, teList, nSpan=10, nChord=10, lower=1.0, scaled=True)
For the volume constraint, the volume is computed by adding up the volumes of the prisms that make up the projected grid as illustrated in the following image (only showing a section for clarity). For the thickness constraints, the distances between the upper and lower projected points are used, as illustrated in the following image. During optimization, these projected points are also moved by the FFD, just like the wing surface, and are used again to calculate the thicknesses and volume for the new designs. More information on the options can be found in the pyGeo docs or by looking at the pyGeo source code.
The LeTe constraints (short for Leading edge/Trailing edge constraints) are linear constraints based on the FFD control points.
When we have both twist and local shape variables, we want to prevent the local shape variables from creating a shearing twist.
This is done by constraining the upper and lower FFD control points on the leading and trailing edges to move in opposite directions.
Note that the LeTe constraint is not related to the leList and teList points discussed above.
# Le/Te constraints
DVCon.addLeTeConstraints(0, "iLow")
DVCon.addLeTeConstraints(0, "iHigh")
if comm.rank == 0:
# Only make one processor do this
DVCon.writeTecplot(os.path.join(args.output, "constraints.dat"))
In this script DVCon.writeTecplot will save a file named constraints.dat which can be opened with Tecplot to visualize and check these constraints.
Since this is added here, before the commands that run the optimization, the file will correspond to the initial geometry.
The following image shows the constraints visualized with the wing surface superimposed.
This command can also be added at the end of the script to visualize the final constraints.
Mesh warping set-up
meshOptions = {"gridFile": args.gridFile}
mesh = USMesh(options=meshOptions, comm=comm)
CFDSolver.setMesh(mesh)
Optimization callback functions
First we must set up a callback function and a sensitivity function for each processor set.
In this case cruiseFuncs and cruiseFuncsSens belong to the cruise processor set.
Then we need to set up an objCon function, which is used to create abstract functions of other functions.
def cruiseFuncs(x):
if comm.rank == 0:
print(x)
# Set design vars
DVGeo.setDesignVars(x)
ap.setDesignVars(x)
# Run CFD
CFDSolver(ap)
# Evaluate functions
funcs = {}
DVCon.evalFunctions(funcs)
CFDSolver.evalFunctions(ap, funcs)
CFDSolver.checkSolutionFailure(ap, funcs)
if comm.rank == 0:
print(funcs)
return funcs
def cruiseFuncsSens(x, funcs):
funcsSens = {}
DVCon.evalFunctionsSens(funcsSens)
CFDSolver.evalFunctionsSens(ap, funcsSens)
CFDSolver.checkAdjointFailure(ap, funcsSens)
return funcsSens
def objCon(funcs, printOK):
# Assemble the objective and any additional constraints:
funcs["obj"] = funcs[ap["cd"]]
funcs["cl_con_" + ap.name] = funcs[ap["cl"]] - 0.5
if printOK:
print("funcs in obj:", funcs)
return funcs
Now we will explain each of these callback functions.
cruiseFuncs
The input to cruiseFuncs is the dictionary of design variables.
First, we pass this dictionary to DVGeometry and AeroProblem to set their respective design variables.
Then we solve the flow solution given by the AeroProblem with ADflow.
Finally, we fill the funcs dictionary with the function values computed by DVConstraints and ADflow.
The call checkSolutionFailure checks ADflow to see if there was a failure in the solution (could be due to negative volumes or something more sinister).
If there was a failure it changes the fail flag in funcs to True.
The funcs dictionary is the required return.
cruiseFuncsSens
The inputs to cruiseFuncsSens are the design variable and function dictionaries.
Inside cruiseFuncsSens we populate the funcsSens dictionary with the derivatives of each of the functions in cruiseFuncs with respect to all of its dependence variables.
objCon
The main input to the objCon callback function is the dictionary of functions (which is a compilation of all the funcs dictionaries from each of the design points).
Inside objCon, the user can define functionals (or functions of other functions).
For instance, to maximize L/D, you could define the objective function as:
funcs['obj'] = funcs['cl'] / funcs['cd']
The objCon function is processed within the multipoint module and the partial derivatives of any functionals with respect to the input functions are automatically computed with the complex-step method.
This means that the user doesn’t have to worry about computing analytic derivatives for the simple functionals defined in objCon.
The printOK input is a boolean that is False when the complex-step is in process.
Optimization problem
Setting up the optimization problem follows the same format as before, only now we incorporate multiPointSparse.
When creating the instance of the Optimization problem, MP.obj is given as the objective function.
multiPointSparse will take care of calling both cruiseFuncs and objCon to provide the full funcs dictionary to pyOptSparse.
Both AeroProblem and DVGeometry have built-in functions to add all of their respective design variables to the optimization problem.
DVConstraints also has a built-in function to add all constraints to the optimization problem.
The user must manually add any constraints that were defined in objCon.
Finally, we need to tell multiPointSparse which callback functions belong to which processor set.
We also need to provide it with the objCon and the optProb.
The call optProb.printSparsity() prints out the constraint Jacobian at the beginning of the optimization.
# Create optimization problem
optProb = Optimization("opt", MP.obj, comm=comm)
# Add objective
optProb.addObj("obj", scale=1e2)
# Add variables from the AeroProblem
ap.addVariablesPyOpt(optProb)
# Add DVGeo variables
DVGeo.addVariablesPyOpt(optProb)
# Add constraints
DVCon.addConstraintsPyOpt(optProb)
optProb.addCon("cl_con_" + ap.name, lower=0.0, scale=10.0)
# The MP object needs the 'obj' and 'sens' function for each proc set,
# the optimization problem and what the objcon function is:
MP.setProcSetObjFunc("cruise", cruiseFuncs)
MP.setProcSetSensFunc("cruise", cruiseFuncsSens)
MP.setObjCon(objCon)
MP.setOptProb(optProb)
optProb.printSparsity()
Run optimization
To finish up, we choose the optimizer and then run the optimization.
# Set up optimizer
if args.opt == "SLSQP":
optOptions = {"IFILE": os.path.join(args.output, "SLSQP.out")}
elif args.opt == "SNOPT":
optOptions = {
"Major feasibility tolerance": 1e-4,
"Major optimality tolerance": 1e-4,
"Hessian full memory": None,
"Function precision": 1e-8,
"Print file": os.path.join(args.output, "SNOPT_print.out"),
"Summary file": os.path.join(args.output, "SNOPT_summary.out"),
"Major iterations limit": 1000,
}
elif args.opt == "IPOPT":
optOptions = {
"limited_memory_max_history": 1000,
"print_level": 5,
"tol": 1e-6,
"acceptable_tol": 1e-5,
"max_iter": 300,
"start_with_resto": "yes",
}
optOptions.update(args.optOptions)
opt = OPT(args.opt, options=optOptions)
# Run Optimization
sol = opt(optProb, MP.sens, storeHistory=os.path.join(args.output, "opt.hst"))
if comm.rank == 0:
print(sol)
Note
The complete set of options for SNOPT can be found in the pyOptSparse documentation.
It is useful to remember that you can include major iterations information in the history file by providing the proper options.
It has been observed that the _print and _summary files occasionally fail to be updated, possibly due to unknown hardware issues on GreatLakes.
The problem is not common, but if you want to avoid losing this information, you might back it up in the history file.
This would allow you monitor the optimization even if the _print and _summary files are not being updated.
Note that the size of the history file will increase due to this additional data.
Run it yourself!
To run the script, use the mpirun and place the total number of processors after the -np argument
mpirun -np 4 python aero_opt.py
You can follow the progress of the optimization using OptView, as explained in pyOptSparse.
Terminal output
An important step of verifyting the optimization setup is to check the sparsity structure of the constraint Jacobian:
+-------------------------------------------------------------------------------+
| Sparsity structure of constraint Jacobian |
+-------------------------------------------------------------------------------+
alpha_wing (1) twist (7) local (96)
+---------------------------------------+
DVCon1_volume_constraint_0 (1) | | X | X |
-----------------------------------------
DVCon1_thickness_constraints_0 (100) | | X | X |
-----------------------------------------
DVCon1_lete_constraint_0_local(L) (8) | | | X |
-----------------------------------------
DVCon1_lete_constraint_1_local(L) (8) | | | X |
-----------------------------------------
cl_con_wing (1) | X | X | X |
+---------------------------------------+
Postprocessing the solution output
All output is found in the output directory.
The naming scheme of the files follows in general <name>_<iter>_<type>.<ext>, where the <name> is aeroproblem, <iter> is the function evaluation number, <type> is the solution type (surface, volume, lift, slices, ect.), and <ext> is the file extension.
The solution files (.dat, .cgns or .plt) can be viewed in the Tecplot.
A contour plot of the pressure coefficient compared with the surface solution from the Analysis with ADflow is shown below.
Similarly, as done in Analysis with ADflow, the lift and slice files (.dat) are used to compare the spanwise normalized lift, compared to elliptical lift, the twist distribution, and t/c.
For the slice file, here the normalized airfoil shape and pressure coefficient are shown.
The optimized design achieves an elliptical lift distribution, and shows a more even and gradual pressure distribution as shown by the contour and section plots.
Further, the optimized twist distribution, as expected, demonstrates higher lift in the outboard section of the wing.
Finally the optimization history can be viewed either by parsing the database or using OptView. Here the former is done showing the major iterations of the history, objective, and some of the design variables and constraints. To produce the figure, run the accompanying postprocessing script as shown below.
python plot_optHist.py --histFile output/opt.hst --outputDir output/
