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:

minimize
\(C_D\)
with respect to
10 shape variables
subject to
\(C_L = 0.5\)

The shape variables are controlled by the FFD points specified in the FFD file.

Files

Copy the FFD file, ffd.xyz, and the CGNS mesh file, n0012.cgns, generated previously, into the directory.

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

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", "cmz", "yplus"],
    "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)
CFDSolver.addLiftDistribution(200, "z")

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:

nCycles

If the analysis doesn’t converge, this can be increased.

nkswitchtol

If 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.

NKSubSpaceSize

Decreasing this parameter will decrease memory usage when in the NK range. Only change this if there are memory issues when dealing with larger meshes.

writeSurfaceSolution and writeVolumeSolution

If you want to view the surface or volume solutions at the end of each analysis, these parameters can be set to True.

We also add a single lift distribution with 200 sampling points using addLiftDistribution, meaning that these values are written to a text file after every iteration. Similarly, addSlices writes the airfoil coordinates of the specified slice 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.

# Create DVGeometry object
FFDFile = "ffd.xyz"

DVGeo = DVGeometry(FFDFile)
DVGeo.addGeoDVLocal("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

Note: This section is also the same as the corresponding section in aircraft optimization.

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.

The volume and thickness constraints are set up by creating a 2D grid of points which is projected through the planform of the wing. For the volume constraint, the 2D grid is transformed into a 3D grid bounded by the surface of the wing. The volume is computed by adding up the volumes of the cells that make up the 3D grid. For the thickness constraints, the nodes of the 2D grid are projected to the upper and lower surface of the wing. The thickness for a given node is the difference between its upper and lower projections.

The LeTe constraints (short for Leading edge/Trailing edge) 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 that the upper and lower nodes on the leading and trailing edges must move in opposite directions.


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=0.064837137176294343, upper=0.07, scaled=False)
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)

The parameters lower and upper are thickness or volume bounds relative to the original airfoil. For example, the lower bound of 0.1 in the thickness constraint means that the optimized airfoil can be no thinner than 10% of the original airfoil.

Mesh warping set-up

meshOptions = {"gridFile": args.gridFile}

mesh = USMesh(options=meshOptions, comm=comm)
CFDSolver.setMesh(mesh)

Optimization callback functions

Note: This section is also the same as the corresponding section in aircraft optimization.

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. This should be similar for all single-point optimizations.

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


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).

Optimization problem

Note: This section is also the same as the corresponding section in aircraft optimization.

To set up the optimization problem, 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=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,
        "Difference interval": 1e-3,
        "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:

$ mkdir output
$ 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.