Single Point Optimization
Introduction
We will now proceed to optimizing a NACA0012 airfoil for a given set of flow conditions. It is very similar to a wing optimization. The optimization problem is defined as:
The shape variables are controlled by the FFD points specified in the FFD file.
Files
Navigate to the directory airfoilopt/singlepoint in your tutorial folder.
Copy the FFD file, ffd.xyz, and the CGNS mesh file, n0012.cgns, generated previously, into the directory:
cp ../mesh/n0012.cgns .
cp ../ffd/ffd.xyz .
Create the following empty runscript in the current directory:
airfoil_opt.py
Dissecting the aerodynamic optimization script
Open the file airfoil_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 numpy as np
import argparse
import ast
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
Adding command line arguments
This is a convenience feature that allows the user to pass in command line arguments to the script. Four options are provided:
Output directory
Optimizer to be used
Grid file to be used
Optimizer options
# 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=["SLSQP", "SNOPT"])
parser.add_argument("--gridFile", type=str, default="n0012.cgns")
parser.add_argument("--optOptions", type=ast.literal_eval, default={}, help="additional optimizer options to be added")
args = parser.parse_args()
Specifying parameters for the optimization
Several conditions for the optimization are specified at the beginning of the script. These include the coefficient of lift constraint value, Mach number, and altitude to indicate flow conditions.
# cL constraint
mycl = 0.5
# angle of attack
alpha = 1.5
# mach number
mach = 0.75
# cruising altitude
alt = 10000
The angle of attack serves as the initial value for the optimization and should not affect the optimized result.
Creating processor sets
Allocating sets of processors for different analyses can be helpful for multiple design points, but this is a single-point optimization, so only one point is added.
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)
ADflow set-up
The ADflow set-up looks similar to the aerodynamic analysis script.
aeroOptions = {
# Common Parameters
"gridFile": args.gridFile,
"outputDirectory": args.output,
# Physics Parameters
"equationType": "RANS",
"smoother": "DADI",
"MGCycle": "sg",
"nCycles": 20000,
"monitorvariables": ["resrho", "cl", "cd"],
"useNKSolver": True,
"useanksolver": True,
"nsubiterturb": 10,
"liftIndex": 2,
"infchangecorrection": True,
# Convergence Parameters
"L2Convergence": 1e-15,
"L2ConvergenceCoarse": 1e-4,
# Adjoint Parameters
"adjointSolver": "GMRES",
"adjointL2Convergence": 1e-12,
"ADPC": True,
"adjointMaxIter": 5000,
"adjointSubspaceSize": 400,
"ILUFill": 3,
"ASMOverlap": 3,
"outerPreconIts": 3,
"NKSubSpaceSize": 400,
"NKASMOverlap": 4,
"NKPCILUFill": 4,
"NKJacobianLag": 5,
"nkswitchtol": 1e-6,
"nkouterpreconits": 3,
"NKInnerPreConIts": 3,
"writeSurfaceSolution": False,
"writeVolumeSolution": False,
"frozenTurbulence": False,
"restartADjoint": False,
}
# Create solver
CFDSolver = ADFLOW(options=aeroOptions, comm=comm)
As it is, the options specified above allow for a good convergence of NACA0012 airfoil analysis, but may not converge for other airfoils. Some useful options to adjust are:
nCyclesIf the analysis doesn’t converge, this can be increased.
nkswitchtolIf the analysis stops converging during NK (Newton-Krylov), this might mean that it is still outside of the radius of convergence of the NK method. The parameter should then be lowered.
NKSubSpaceSizeDecreasing this parameter will decrease memory usage when in the NK range. Only change this if there are memory issues when dealing with larger meshes.
writeSurfaceSolutionandwriteVolumeSolutionIf you want to view the surface or volume solutions at the end of each analysis, these parameters can be set to True.
We also use addSlices to write the airfoil coordinates and \(c_p\) distribution to a text file.
Set the AeroProblem
We add angle of attack as a design variable and set up the AeroProblem using given flow conditions.
ap = AeroProblem(name="fc", alpha=alpha, mach=mach, altitude=alt, areaRef=1.0, chordRef=1.0, evalFuncs=["cl", "cd"])
# Add angle of attack variable
ap.addDV("alpha", value=alpha, lower=0, upper=10.0, scale=1.0)
Geometric parametrization
The set-up for DVGeometry is simpler for an airfoil since it doesn’t involve span-wise variables such as twist, dihedral, or taper. As a result, we also don’t need to set up a reference axis. The only DVs we have are local shape variables that control the vertical movements of each individual FFD node. Note that since we have to work with a 3D problem, this in fact has twice as many DVs as we’d like—the shape of the two airfoil sections should remain the same. This is addressed by adding linear constraints in the following section.
# Create DVGeometry object
FFDFile = "ffd.xyz"
DVGeo = DVGeometry(FFDFile)
DVGeo.addLocalDV("shape", lower=-0.05, upper=0.05, axis="y", scale=1.0)
span = 1.0
pos = np.array([0.5]) * span
CFDSolver.addSlices("z", pos, sliceType="absolute")
# Add DVGeo object to CFD solver
CFDSolver.setDVGeo(DVGeo)
The local design variable shape is added.
Geometric constraints
This section is very similar to the corresponding section for the wing optimization.
The only difference is that, we must add a set of linear constraints such that the shape deformations on one side of the airfoil mirrors that of the other.
This is accomplished with a call to addLinearConstraintsShape.
DVCon = DVConstraints()
DVCon.setDVGeo(DVGeo)
# Only ADflow has the getTriangulatedSurface Function
DVCon.setSurface(CFDSolver.getTriangulatedMeshSurface())
# Le/Te constraints
lIndex = DVGeo.getLocalIndex(0)
indSetA = []
indSetB = []
for k in range(0, 1):
indSetA.append(lIndex[0, 0, k]) # all DV for upper and lower should be same but different sign
indSetB.append(lIndex[0, 1, k])
for k in range(0, 1):
indSetA.append(lIndex[-1, 0, k])
indSetB.append(lIndex[-1, 1, k])
DVCon.addLeTeConstraints(0, indSetA=indSetA, indSetB=indSetB)
# DV should be same along spanwise
lIndex = DVGeo.getLocalIndex(0)
indSetA = []
indSetB = []
for i in range(lIndex.shape[0]):
indSetA.append(lIndex[i, 0, 0])
indSetB.append(lIndex[i, 0, 1])
for i in range(lIndex.shape[0]):
indSetA.append(lIndex[i, 1, 0])
indSetB.append(lIndex[i, 1, 1])
DVCon.addLinearConstraintsShape(indSetA, indSetB, factorA=1.0, factorB=-1.0, lower=0, upper=0)
le = 0.0001
leList = [[le, 0, le], [le, 0, 1.0 - le]]
teList = [[1.0 - le, 0, le], [1.0 - le, 0, 1.0 - le]]
DVCon.addVolumeConstraint(leList, teList, 2, 100, lower=1, scaled=True)
DVCon.addThicknessConstraints2D(leList, teList, 2, 100, lower=0.1, upper=3.0)
if comm.rank == 0:
fileName = os.path.join(args.output, "constraints.dat")
DVCon.writeTecplot(fileName)
Mesh warping set-up
meshOptions = {"gridFile": args.gridFile}
mesh = USMesh(options=meshOptions, comm=comm)
CFDSolver.setMesh(mesh)
Optimization callback functions
This section is also the same as the corresponding section in aircraft optimization.
def cruiseFuncs(x):
if MPI.COMM_WORLD.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 MPI.COMM_WORLD.rank == 0:
print(funcs)
return funcs
def cruiseFuncsSens(x, funcs):
funcsSens = {}
DVCon.evalFunctionsSens(funcsSens)
CFDSolver.evalFunctionsSens(ap, funcsSens)
CFDSolver.checkAdjointFailure(ap, funcsSens)
if MPI.COMM_WORLD.rank == 0:
print(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"]] - mycl
if printOK:
print("funcs in obj:", funcs)
return funcs
Optimization problem
This section is also the same as the corresponding section in aircraft optimization.
# Create optimization problem
optProb = Optimization("opt", MP.obj, comm=MPI.COMM_WORLD)
# Add objective
optProb.addObj("obj", scale=1e4)
# 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, upper=0.0, scale=1.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
# 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"),
}
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 MPI.COMM_WORLD.rank == 0:
print(sol)
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 airfoil_opt.py | tee output.txt
The command tee saves the text outputs of the optimization to the specified text file.
You can follow the progress of the optimization using OptView, as explained in pyOptSparse.
