# 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*

*with respect to*

*subject to*

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.