"""
CMPLXFOIL
CMPLXFOIL is a wrapper for Mark Drela's XFOIL code. The purpose of this
class is to provide an easy to use wrapper for xfoil for intergration
into other projects. Both real and complex versions of xfoil can be used.
The version uses the MDO Lab's BaseSolver baseclass so it can be used
within the MACH-Aero optimization framework.
Developers:
-----------
- Eytan Adler (EA)
- Alasdair Gray (AG)
History
-------
v. 1.0 - Initial Class Creation (EA + AG, 2022)
"""
# =============================================================================
# Standard Python modules
# =============================================================================
import os
import time
from copy import copy, deepcopy
import pickle as pkl
from collections.abc import Iterable
import warnings
# =============================================================================
# Other Python modules
# =============================================================================
import numpy as np
from baseclasses import BaseSolver
from . import MExt
# Handle plotting imports
plt = None
try:
from matplotlib import pyplot
import matplotlib.lines as mpl_lines
plt = pyplot
except ImportError:
pass
# If the matplotlib import worked, try niceplots
if plt:
try:
import niceplots
plt.style.use(niceplots.get_style())
colors = niceplots.get_colors()
color = colors["Blue"]
cpUpColor = colors["Blue"]
cpLowColor = colors["Red"]
except ImportError:
color = "b"
cpUpColor = "b"
cpLowColor = "r"
[docs]
class CMPLXFOIL(BaseSolver):
"""
CMPLXFOIL Class Initialization
Parameters
----------
fileName : str
Filename of DAT file to read in
options : dict of option-value pairs, optional
Options for the solver. Available options can be found
in the Options section of the documentation or the
options.yaml file in the docs directory.
debug : bool, optional
Set this flag to true when debugging with a symbolic
debugger. The MExt module deletes the copied .so file when not
required which causes issues debugging, by default False
"""
def __init__(self, fileName, options={}, debug=False):
# Load the compiled module using MExt, allowing multiple imports
curDir = os.path.basename(os.path.dirname(os.path.realpath(__file__)))
time.sleep(0.1) # BOTH of these sleeps are necessary for some reason!
self.xfoil = MExt.MExt("libcmplxfoil", curDir, debug=debug)._module
time.sleep(0.1) # BOTH of these sleeps are necessary for some reason!
self.xfoil_cs = MExt.MExt("libcmplxfoil_cs", curDir, debug=debug)._module
super().__init__("CMPLXFOIL", "Panel Method", defaultOptions=self._getDefaultOptions(), options=options)
if self.getOption("xTrip").shape != (2,):
raise ValueError("xTrip option shape must be (2,)")
# Read in DAT file and create coordinates. DVGeo needs 3D coordinates
# so keep the third dimension as dummy coordinates.
self.DATFileName = fileName
self.coords = self._readDat(fileName)
self.coords = np.hstack((self.coords, np.zeros((self.coords.shape[0], 1))))
self.coords0 = self.coords.copy() # initial coordinates (never changes)
self.setCoordinates(self.coords0) # set the initial coordinates
self.setCoordinatesComplex(self.coords0) # set the initial complex coordinates
self.pointSetKwargs = {}
self.curAP = None
self.DVGeo = None
# Dictionary with dictionary of functions for each aero problem
self.funcs = {}
self.funcsComplex = {}
self.functionList = ["cl", "cd", "cm"] # available functions
# When the XFOIL solver is called, slice data is saved (key is the current
# AeroProblem name). In the associated value is a dictionary containing
# "cp_visc_upper": viscous CP on the airfoil's upper surface
# "cp_invisc_upper": inviscid CP on the airfoil's upper surface
# "x_upper": x coordinates of the upper surface CP data
# "y_upper": y coordinates of the upper surface CP data
# "cp_visc_lower": viscous CP on the airfoil's lower surface
# "cp_invisc_lower": inviscid CP on the airfoil's lower surface
# "x_lower": x coordinates of the lower surface CP data
# "y_lower": y coordinates of the lower surface CP data
# "cf_upper": skin friction coefficient on the upper surface
# "x_cf_upper": x coordinates of upper surface skin friction coefficient
# "y_cf_upper": y coordinates of upper surface skin friction coefficient
# "cf_lower": skin friction coefficient on the lower surface
# "x_cf_lower": x coordinates of lower surface skin friction coefficient
# "y_cf_lower": y coordinates of lower surface skin friction coefficient
self.sliceData = {}
self.sliceDataComplex = {}
# Possible AeroProblem design variables (only alpha for CMPLXFOIL)
self.possibleAeroDVs = ["alpha"]
# Figure and axes used by self.plotAirfoil and self._updateAirfoilPlot
self.airfoilFig = None
self.airfoilAxs = None
self.CPlim = None # y limits on the CP plot
self.CFlim = None # y limits on the CF plot
[docs]
def setDVGeo(self, DVGeo, pointSetKwargs=None):
"""
Set the DVGeometry object that will manipulate 'geometry' in
this object. Note that CMPLXFOIL does not **strictly** need a
DVGeometry object, but if optimization with geometric
changes is desired, then it is required.
Parameters
----------
DVGeo : A DVGeometry object.
Object responsible for manipulating the constraints that
this object is responsible for.
pointSetKwargs : dict
Keyword arguments to be passed to the DVGeo addPointSet call.
Useful for DVGeometryMulti, specifying FFD projection tolerances, etc
"""
self.DVGeo = DVGeo
# Get the number of geometry variables
self.numGeoDV = self.DVGeo.getNDV()
# Save kwargs for addPointSet
if pointSetKwargs is not None:
self.pointSetKwargs = pointSetKwargs
[docs]
def setAeroProblem(self, aeroProblem):
"""
Sets the aeroProblem to by used by CMPLXFOIL.
Parameters
----------
aeroProblem : :class:`AeroProblem <baseclasses.problems.pyAero_problem.AeroProblem>` instance
Aero problem to set (gives flight conditions)
"""
ptSetName = "cmplxfoil_%s_coords" % aeroProblem.name
# Tell the user if we are switching aeroProblems
if self.curAP != aeroProblem:
self.pp("+" + "-" * 70 + "+")
self.pp("| Switching to Aero Problem: %-41s|" % aeroProblem.name)
self.pp("+" + "-" * 70 + "+")
# Now check if we have an DVGeo object to deal with:
if self.DVGeo is not None:
# DVGeo appeared and we have not embedded points!
if ptSetName not in self.DVGeo.points:
self.DVGeo.addPointSet(self.coords0, ptSetName, **self.pointSetKwargs)
aeroProblem.ptSetName = ptSetName
# Check if our point-set is up to date:
if not self.DVGeo.pointSetUpToDate(ptSetName):
coords = self.DVGeo.update(ptSetName, config=aeroProblem.name)
self.setCoordinates(coords)
# Initialize the callCounter if it's not already an attribute
try:
aeroProblem.callCounter
except AttributeError:
aeroProblem.callCounter = 0
self.curAP = aeroProblem
[docs]
def setCoordinates(self, coords):
"""
Update the airfoil coordinates and associated point sets.
Parameters
----------
coords : ndarray
New airfoil coordinates (either 2 or 3 columns)
"""
self.coords[:, :2] = coords[:, :2].copy()
self._setCoordinates(self.coords[:, 0], self.coords[:, 1])
[docs]
def getCoordinates(self):
"""
Return the current airfoil coordinates
Returns
-------
coords : ndarray
Airfoil coordinates with each column being (x, y, z)
where z is a dummy value.
"""
return self.coords.copy()
[docs]
def setCoordinatesComplex(self, coords):
"""
Update the complex airfoil coordinates and associated point sets.
Parameters
----------
coords : ndarray
New airfoil coordinates (either 2 or 3 columns)
"""
self._setCoordinates(coords[:, 0].astype(complex), coords[:, 1].astype(complex), setComplex=True)
def __call__(self, aeroProblem, useComplex=False, deriv=False):
"""
Evaluate XFOIL with the current coordinates and flight conditions (from aeroProblem).
Parameters
----------
aeroProblem : :class:`AeroProblem <baseclasses.problems.pyAero_problem.AeroProblem>` instance
Aero problem to set (gives flight conditions)
useComplex : bool
Run XFOIL in complex mode
deriv : bool
Set to True if the call is associated with the derivative calculation. callCounter will be incremented and
writeSolution called only if this is False
"""
if useComplex:
dtype = complex
xfoil = self.xfoil_cs
funcs = self.funcsComplex
sliceData = self.sliceDataComplex
else:
dtype = float
xfoil = self.xfoil
funcs = self.funcs
sliceData = self.sliceData
self.setAeroProblem(aeroProblem)
# Set flight condition and options
xfoil.cr15.reinf1 = aeroProblem.re # Reynolds number
xfoil.cr09.minf1 = aeroProblem.mach # Mach Number set
xfoil.cr09.adeg = aeroProblem.alpha
xfoil.ci04.itmax = self.getOption("maxIters") # Iterations Limit Set
if not np.any(np.isnan(self.getOption("xTrip"))): # NaN is default to not set, otherwise set it
xfoil.cr15.xstrip = self.getOption("xTrip")
# Call XFOIL
xfoil.oper()
# Store results in dictionary for current aero problem
funcs[aeroProblem.name] = {
"cl": dtype(xfoil.cr09.cl),
"cd": dtype(xfoil.cr09.cd),
"cm": dtype(xfoil.cr09.cm),
}
# Pull out and process the pressure and skin friction coefficient data
cpv = xfoil.cr04.cpv # viscous
cpi = xfoil.cr04.cpi # inviscid
x = xfoil.cr05.x
y = xfoil.cr05.y
end_foil_idx = np.argmax(x > self.coords[0, 0]) + 1 # XFOIL includes wake panels, which we don't want
idx_lower_start = end_foil_idx // 2 # first half of data is upper surface
idxUpper = np.argwhere(xfoil.cr15.tau[:, 0] != 0).flatten()
idxLower = np.argwhere(xfoil.cr15.tau[:, 1] != 0).flatten()
iPanUpper = xfoil.ci05.ipan[idxUpper, 0] - 1 # FORTRAN uses 1-based indexing, need to adjust
iPanLower = xfoil.ci05.ipan[idxLower, 1] - 1 # FORTRAN uses 1-based indexing, need to adjust
tauUpper = xfoil.cr15.tau[idxUpper, 0].copy().astype(dtype)
tauLower = xfoil.cr15.tau[idxLower, 1].copy().astype(dtype)
uInf = xfoil.cr09.qinf
xCfUpper = x[iPanUpper].copy().astype(dtype)
yCfUpper = y[iPanUpper].copy().astype(dtype)
xCfLower = x[iPanLower].copy().astype(dtype)
yCfLower = y[iPanLower].copy().astype(dtype)
sliceData[aeroProblem.name] = {
"cp_visc_upper": cpv[:idx_lower_start].copy().astype(dtype),
"cp_invisc_upper": cpi[:idx_lower_start].copy().astype(dtype),
"x_upper": x[:idx_lower_start].copy().astype(dtype),
"y_upper": y[:idx_lower_start].copy().astype(dtype),
"cp_visc_lower": cpv[idx_lower_start:end_foil_idx].copy().astype(dtype),
"cp_invisc_lower": cpi[idx_lower_start:end_foil_idx].copy().astype(dtype),
"x_lower": x[idx_lower_start:end_foil_idx].copy().astype(dtype),
"y_lower": y[idx_lower_start:end_foil_idx].copy().astype(dtype),
"cf_upper": tauUpper / (0.5 * uInf**2),
"x_cf_upper": xCfUpper,
"y_cf_upper": yCfUpper,
"cf_lower": tauLower / (0.5 * uInf**2),
"x_cf_lower": xCfLower,
"y_cf_lower": yCfLower,
}
# Check for failure
self.curAP.solveFailed = self.curAP.fatalFail = xfoil.cl01.lexitflag != 0 or xfoil.cl01.lvconv == 0
# If not a derivative call, increment callCounter
if not deriv:
self.curAP.callCounter += 1
# Write solution files if desired
if not deriv and self.getOption("writeSolution"):
self.writeSolution()
[docs]
def solveCL(
self,
aeroProblem,
CLStar,
alpha0=None,
alphaBound=None,
delta=0.5,
tol=1e-3,
CLalphaGuess=None,
maxIter=20,
useNewton=False,
):
"""Find the angle of attack that gives a target lift coefficient.
Parameters
----------
aeroProblem : pyAero_problem class
The aerodynamic problem to solve
CLStar : float
The desired CL
alpha0 : float, optional
Initial guess for secant search (deg). If None, use the value in the aeroProblem, by default None
alphaBound : float, tuple, list, optional
Bounds for angle of attack, if scalar then value is treated as a +- bound, by default None, in which case
limit is +/-15 deg
delta : float, optional
Initial step direction for secant search, by default 0.5
tol : float, optional
Desired tolerance for CL, by default 1e-3
CLalphaGuess : float, optional
The user can provide an estimate for the lift curve slope in order to accelerate convergence. If the user
supplies a value to this option, it will not use the delta value anymore to select the angle of attack of
the second run. The value should be in 1/deg., by default None
maxIter : int, optional
Maximum number of iterations, by default 20
useNewton : bool, optional
If True, Newton's method will be used where the dCL/dAlpha is computed using complex-step, otherwise the
secant method is used, by default False
Returns
-------
None, but the correct alpha is stored in the aeroProblem
"""
self.setAeroProblem(aeroProblem)
if alpha0 is not None:
aeroProblem.alpha = alpha0
else:
alpha0 = aeroProblem.alpha
if alphaBound is None:
alphaBound = (-15, 15)
elif isinstance(alphaBound, (int, float)):
alphaBound = (-alphaBound, alphaBound)
elif isinstance(alphaBound, Iterable):
alphaBound = (alphaBound[0], alphaBound[1])
else:
raise ValueError(
f'Supplied alphaBound value "{alphaBound}" is not the correct type, must be a scalar or array-like value.'
)
dCLdAlpha = CLalphaGuess
resPrev = None
alphaPrev = None
hasConverged = False
maxRes = 1e2
for _ in range(maxIter):
# Call the solver and compute the "residual"
self.__call__(aeroProblem, deriv=True)
failCheck = {}
self.checkSolutionFailure(aeroProblem, failCheck)
cl = float(self.funcs[aeroProblem.name]["cl"])
print(f"Alpha: {aeroProblem.alpha:>6.3f}, CL: {cl:>7.6f}")
res = cl - CLStar
# Check for convergence and divergence
hasConverged = abs(res) < tol and not failCheck["fail"]
if hasConverged:
print("Converged!")
break
elif abs(res > maxRes):
warnings.warn(f"solveCL diverged (CL = {cl:.2f})", stacklevel=2)
break
# If not converged or diverged, compute the next alpha
if not failCheck["fail"]:
if useNewton:
aeroProblem.alpha += 1e-200 * 1j
self.__call__(aeroProblem, useComplex=True, deriv=True)
dCLdAlpha = np.imag(self.funcsComplex[aeroProblem.name]["cl"]) * 1e200
aeroProblem.alpha = np.real(aeroProblem.alpha)
else:
# If using secant, we can only compute dCLdAlpha from the second iteration onwards
if resPrev is not None:
dCLdAlpha = (res - resPrev) / (aeroProblem.alpha - alphaPrev)
resPrev = res
alphaPrev = aeroProblem.alpha
# Update the alpha either using dCLdAlpha or delta, or backtracking if something went wrong
if dCLdAlpha == 0.0 or failCheck["fail"]:
aeroProblem.alpha *= 0.9
elif dCLdAlpha is not None:
aeroProblem.alpha = np.clip(aeroProblem.alpha - res / dCLdAlpha, alphaBound[0], alphaBound[1])
else:
aeroProblem.alpha = aeroProblem.alpha + delta
if not hasConverged:
print("Did not converge :(")
return hasConverged
[docs]
def writeSolution(self, outputDir=None, baseName=None, number=None):
"""This is a generic shell function that potentially writes the various output files. The intent is that the
user or calling program can call this file and CMPLXFOIL write all the files that the user has defined. It is
recommended that this function is used along with the associated logical flags in the options to determine the
desired writing procedure.
Parameters
----------
outputDir : str
Use the supplied output directory
baseName : str
Use this supplied string for the base filename. Typically only used from an external solver.
number : int
Use the user supplied number to index solution. Again, only typically used from an external solver.
"""
if outputDir is None:
outputDir = self.getOption("outputDirectory")
if baseName is None:
baseName = self.curAP.name
# Add a number to the filename, either from the user or from the current callCounter
if number is not None:
baseName = baseName + "_%3.3d" % number
else:
if self.getOption("numberSolutions"):
baseName = baseName + "_%3.3d" % self.curAP.callCounter
# Join to get the actual filename root
base = os.path.join(outputDir, baseName)
if self.getOption("writeCoordinates"):
self.writeCoordinates(base)
if self.getOption("writeSliceFile"):
self.writeSlice(base)
if self.getOption("plotAirfoil"):
self.plotAirfoil()
[docs]
def writeCoordinates(self, fileName):
"""Write dat file with the current coordinates.
Parameters
----------
fileName : str
File name for saved dat file (".dat" will be automatically appended).
"""
fileName += ".dat"
with open(fileName, "w") as f:
for i in range(self.coords.shape[0]):
f.write(str(round(self.coords[i, 0], 12)) + "\t\t" + str(round(self.coords[i, 1], 12)) + "\n")
[docs]
def writeSlice(self, fileName):
"""Write pickle file containing the sliceData dictionary. The data can be
accessed using the AeroProblem name as the key. Within that is a dictionary containing
* Pressure coefficient data on the upper surface
* ``"cp_visc_upper"``: viscous CP on the airfoil's upper surface
* ``"cp_invisc_upper"``: inviscid CP on the airfoil's upper surface
* ``"x_upper"``: x coordinates of the upper surface CP data
* ``"y_upper"``: y coordinates of the upper surface CP data
* Pressure coefficient data on the lower surface
* ``"cp_visc_lower"``: viscous CP on the airfoil's lower surface
* ``"cp_invisc_lower"``: inviscid CP on the airfoil's lower surface
* ``"x_lower"``: x coordinates of the lower surface CP data
* ``"y_lower"``: y coordinates of the lower surface CP data
* Skin friction coefficient data on the upper surface
* ``"cf_upper"``: skin friction coefficient on the upper surface
* ``"x_cf_upper"``: x coordinates of upper surface skin friction coefficient
* ``"y_cf_upper"``: y coordinates of upper surface skin friction coefficient
* Skin friction coefficient data on the lower surface
* ``"cf_lower"``: skin friction coefficient on the lower surface
* ``"x_cf_lower"``: x coordinates of lower surface skin friction coefficient
* ``"y_cf_lower"``: y coordinates of lower surface skin friction coefficient
Parameters
----------
fileName : str
File name for saved pkl file (".pkl" will be automatically appended).
"""
fileName += ".pkl"
with open(fileName, "wb") as f:
pkl.dump(self.sliceData, f, protocol=pkl.HIGHEST_PROTOCOL)
[docs]
def checkSolutionFailure(self, aeroProblem, funcs):
"""Take in a an aeroProblem and check for failure.
Then append the fail flag in funcs. Information regarding whether or not the last analysis with the aeroProblem
was sucessful is included. This information is included as "funcs['fail']". If the 'fail' entry already exits in
the dictionary the following operation is performed:
funcs['fail'] = funcs['fail'] or <did this problem fail>
In other words, if any one problem fails, the funcs['fail'] entry will be True. This information can then be
used directly in multiPointSparse. For direct interface with pyOptSparse the fail flag needs to be returned
separately from the funcs.
Parameters
----------
aeroProblem : pyAero_problem class
The aerodynamic problem to get the solution for
funcs : dict
Dictionary into which the functions are saved.
"""
self.setAeroProblem(aeroProblem)
# We also add the fail flag into the funcs dictionary. If fail is already there, we just logically 'or' what was
# there. Otherwise we add a new entry.
failFlag = self.curAP.solveFailed or self.curAP.fatalFail
if "fail" in funcs:
funcs["fail"] = funcs["fail"] or failFlag
else:
funcs["fail"] = failFlag
[docs]
def checkAdjointFailure(self, aeroProblem, funcsSens):
"""
Pass through to checkSolutionFailure to maintain the same interface as ADflow.
This checks if the primal solve fails and can be called when the sensitivity
is being evaluated (either through FD or CS).
Parameters
----------
aeroProblem : pyAero_problem class
The aerodynamic problem to to get the solution for
funcsSens : dict
Dictionary into which the functions are saved.
"""
self.checkSolutionFailure(aeroProblem, funcsSens)
[docs]
def evalFunctions(self, aeroProblem, funcs, evalFuncs=None, ignoreMissing=False):
"""
This is the main routine for returning useful information from CMPLXFOIL.
The functions corresponding to the strings in ``evalFuncs`` are evaluated
and updated into the provided dictionary.
Parameters
----------
aeroProblem : :class:`AeroProblem <baseclasses.problems.pyAero_problem.AeroProblem>` instance
Aero problem from which to pull evalFuncs and flight conditions.
funcs : dict
Dictionary into which the functions are saved.
evalFuncs : iterable object containing strings
If not none, use these functions to evaluate.
ignoreMissing : bool
Flag to supress checking for a valid function. Please use
this option with caution.
"""
self.setAeroProblem(aeroProblem)
if evalFuncs is None:
evalFuncs = sorted(list(self.curAP.evalFuncs))
else:
evalFuncs = sorted(list(evalFuncs))
# Make the functions lower case
evalFuncs = [s.lower() for s in evalFuncs]
returnFuncs = {}
for f in evalFuncs:
# If it's not in the list of available functions
if f not in self.functionList:
# Either throw an error (if requested) or skip it
if not ignoreMissing:
raise ValueError(f'Supplied function "{f}" is not in the available functions {self.functionList}.')
else:
returnFuncs[aeroProblem.name + "_" + f] = self.funcs[aeroProblem.name][f]
funcs.update(returnFuncs)
# Set the AeroProblem's function names
for name in evalFuncs:
self.curAP.funcNames[name] = self.curAP.name + "_" + name
[docs]
def evalFunctionsSens(self, aeroProblem, funcsSens, evalFuncs=None, mode="CS", h=None):
"""
Evaluate the sensitivity of the desired functions given in
iterable object, 'evalFuncs' and add them to the dictionary
'funcSens'. The keys in the 'funcsSens' dictionary will be have an
``<ap.name>_`` prepended to them.
Parameters
----------
funcsSens : dict
Dictionary into which the function derivatives are saved
evalFuncs : iterable object containing strings
The additional functions the user wants returned that are
not already defined in the aeroProblem
mode : str ["FD" or "CS"]
Specifies how the jacobian vector products will be computed
h : float
Step sized used when the mode is "FD" or "CS" (must be complex
if mode = "CS")
"""
# This is the one and only gateway to the getting derivatives
# out of xfoil. If you want a derivative, it should come from
# here. Thank you.
self.setAeroProblem(aeroProblem)
if evalFuncs is None:
evalFuncs = sorted(list(self.curAP.evalFuncs))
else:
evalFuncs = sorted(list(evalFuncs))
# Make the functions lower case
evalFuncs = [s.lower() for s in evalFuncs]
# Get design variables
DVs = self.DVGeo.getValues()
for dv in self.curAP.DVs.values():
DVs[dv.key] = np.atleast_1d(dv.value)
# Preallocate the funcsSens dictionary with zeros for the desired sensitivities
for f in evalFuncs:
funcsSens[self.curAP.name + "_" + f] = {}
for dvName, dvVal in DVs.items():
if isinstance(dvVal, np.ndarray):
funcsSens[self.curAP.name + "_" + f][dvName] = np.zeros(dvVal.shape, dtype=float)
else:
funcsSens[self.curAP.name + "_" + f][dvName] = np.zeros(1)
# Loop over design variables and compute derivatives of each
for dvName, dvVal in DVs.items():
lenDV = 1 if not isinstance(dvVal, np.ndarray) else len(dvVal)
for i in range(lenDV):
# Compute the seed for the finite difference/complex step
seed = {dvName: np.zeros(dvVal.shape, dtype=float)}
seed[dvName][i] = 1.0
# Compute the sensitivity
sens = self.computeJacobianVectorProductFwd(xDvDot=seed, mode=mode, h=h)
for f in evalFuncs:
funcsSens[self.curAP.name + "_" + f][dvName][i] = sens[f]
# Check that the solution converged
self.checkSolutionFailure(self.curAP, funcsSens)
# Append "_" + (current aero problem name) to any design variables associated with the aero problem
for f in evalFuncs:
func = self.curAP.name + "_" + f
for dvName in DVs.keys():
if dvName in self.possibleAeroDVs and dvName in funcsSens[func]:
funcsSens[func][dvName + "_" + self.curAP.name] = funcsSens[func][dvName]
del funcsSens[func][dvName]
[docs]
def computeJacobianVectorProductFwd(self, xDvDot=None, xSDot=None, mode="CS", h=None):
"""This the main Python gateway for producing forward mode jacobian
vector products. They are computed using either the complex step or finite
difference method. This function is not generally called by the user but
rather internally or from another solver. A DVGeo object must be set
for this routine.
Parameters
----------
xDvDot : dict
Perturbation on the design variables
xSDot : numpy array
Perturbation on the surface
mode : str ["FD" or "CS"]
Specifies how the jacobian vector products will be computed
h : float
Step sized used when the mode is "FD" or "CS" (must be complex
if mode = "CS"), by default 1e-6 for FD and 1e-200j for CS
Returns
-------
dict
Jacobian vector product of evalFuncs given perturbation
"""
if xDvDot is None and xSDot is None:
raise ValueError("xDvDot and xSDot cannot both be None")
if mode not in ["FD", "CS"]:
raise ValueError(f'Jacobian vector product mode "{mode}" invalid. Must be either "FD" or "CS"')
if self.DVGeo is None:
raise ValueError("DVGeo object must be added with setDVGeo before calling computeJacobianVectorProductFwd")
geoDVs = list(self.DVGeo.getValues().keys())
possibleDVs = self.possibleAeroDVs + geoDVs
for DV in xDvDot.keys():
if DV not in possibleDVs:
raise ValueError(f'Perturbed design variable "{DV}" is not valid')
if h is None:
if mode == "FD":
h = 1e-6
elif mode == "CS":
h = 1e-200j
if mode == "FD":
orig_funcs = deepcopy(self.funcs[self.curAP.name])
# Process the Xs perturbation
if xSDot is None:
xsdot = np.zeros_like(self.coords0)
else:
xsdot = xSDot
# For the geometric xDvDot perturbation we accumulate into the
# already existing (and possibly nonzero) xsdot and xvdot
geoxDvDot = {k: xDvDot[k] for k in geoDVs if k in xDvDot}
if xDvDot is not None and self.DVGeo is not None:
xsdot += self.DVGeo.totalSensitivityProd(geoxDvDot, self.curAP.ptSetName, config=self.curAP.name).reshape(
xsdot.shape
)
# Perturb the coordinates
orig_coords = self.getCoordinates()
if mode == "FD":
self.setCoordinates(orig_coords + xsdot * h)
else:
self.setCoordinatesComplex(orig_coords + xsdot * h)
orig_alpha = copy(self.curAP.alpha)
if "alpha" in xDvDot.keys():
self.curAP.alpha += xDvDot["alpha"] * h
self.__call__(self.curAP, useComplex=mode == "CS", deriv=True)
# Compute the Jacobian vector products
jacVecProd = {}
for f in self.functionList:
if mode == "FD":
jacVecProd[f] = (self.funcs[self.curAP.name][f] - orig_funcs[f]) / h
else:
jacVecProd[f] = np.imag(self.funcsComplex[self.curAP.name][f]) / np.imag(h)
# Reset the perturbed variables
self.curAP.alpha = orig_alpha
if mode == "FD":
self.setCoordinates(orig_coords)
self.funcs[self.curAP.name] = orig_funcs
else:
self.setCoordinatesComplex(orig_coords)
return jacVecProd
[docs]
def getTriangulatedMeshSurface(self, offsetDist=1.0):
"""
This function returns a pyGeo surface. The intent is
to use this for DVConstraints.
.. note::
This method requires the pyGeo library
Parameters
----------
offsetDist : float
Distance to extrude airfoil (same units as airfoil coordinates)
Returns
-------
pyGeo surface
Extruded airfoil surface
"""
try:
from pygeo.pyGeo import pyGeo
except ImportError:
raise ImportError("pygeo is required to use getTriangulatedMeshSurface")
airfoil_list = [self.DATFileName] * 2
naf = len(airfoil_list)
# Airfoil leading edge positions
x = [0.0, 0.0]
y = [0.0, 0.0]
z = [0.0, offsetDist]
offset = np.zeros((naf, 2)) # x-y offset applied to airfoil position before scaling
# Airfoil rotations
rot_x = [0.0, 0.0]
rot_y = [0.0, 0.0]
rot_z = [0.0, 0.0]
# Airfoil scaling
scale = [1.0, 1.0] # scaling factor on chord lengths
# Find the trailing edge thickness
teThickness = self.coords[0, 1] - self.coords[-1, 1]
return pyGeo(
"liftingSurface",
xsections=airfoil_list,
scale=scale,
offset=offset,
x=x,
y=y,
z=z,
rotX=rot_x,
rotY=rot_y,
rotZ=rot_z,
bluntTe=True,
squareTeTip=True,
teHeight=teThickness,
)
[docs]
def plotAirfoil(self, fileName=None, showPlot=True):
"""
Plots the current airfoil and returns the figure.
Parameters
----------
fileName : str, optional
FileName to save to, if none specified it will
show the plot with plt.show()
showPlot : bool, optional
Pop open the plot, by default True
Returns
-------
matplotlib figure
Figure with airfoil plotted to it
list of matplotlib axes
List of matplotlib axes for CP and airfoil plots (in that order)
"""
if plt is None:
raise ImportError("matplotlib is required to use the plotting tools")
if self.airfoilFig is None:
# Get data to plot
x = self.coords[:, 0]
y = self.coords[:, 1]
CPUpper = self.sliceData[self.curAP.name]["cp_visc_upper"]
CPUpperInvisc = self.sliceData[self.curAP.name]["cp_invisc_upper"]
xUpper = self.sliceData[self.curAP.name]["x_upper"]
CPLower = self.sliceData[self.curAP.name]["cp_visc_lower"]
CPLowerInvisc = self.sliceData[self.curAP.name]["cp_invisc_lower"]
xLower = self.sliceData[self.curAP.name]["x_lower"]
CFUpper = self.sliceData[self.curAP.name]["cf_upper"]
xCFUpper = self.sliceData[self.curAP.name]["x_cf_upper"]
CFLower = self.sliceData[self.curAP.name]["cf_lower"]
xCFLower = self.sliceData[self.curAP.name]["x_cf_lower"]
# Inverted CP axis
stackedCP = np.hstack((CPUpper, CPUpperInvisc, CPLower, CPLowerInvisc))
stackedCP = stackedCP[np.isfinite(stackedCP)]
stackedCF = np.hstack((CFUpper, CFLower))
stackedCF = stackedCF[np.isfinite(stackedCF)]
self.CPlim = [max(stackedCP) + 0.2, min(stackedCP) - 0.2]
self.CFlim = [min(stackedCF) - 0.002, max(stackedCF) + 0.002]
self.xlimFoil = [min(x) - 0.01, max(x) + 0.01]
self.ylimFoil = [min(y) - 0.01, max(y) + 0.01]
fig, axs = plt.subplots(nrows=3, ncols=1, sharex=True, figsize=[10, 13])
iAxsCP = 0
iAxsCF = 1
iAxsFoil = 2
plt.ion()
if showPlot:
plt.show()
# Plot the CP on the upper axis
axs[iAxsCP].plot(xUpper, CPUpper, color="k", zorder=-1, alpha=0.15)
axs[iAxsCP].plot(xLower, CPLower, color="k", zorder=-1, alpha=0.15)
axs[iAxsCP].plot(xUpper, CPUpper, color=cpUpColor)
axs[iAxsCP].plot(xLower, CPLower, color=cpLowColor)
axs[iAxsCP].plot(xUpper, CPUpperInvisc, "--", color=cpUpColor, linewidth=1.0)
axs[iAxsCP].plot(xLower, CPLowerInvisc, "--", color=cpLowColor, linewidth=1.0)
axs[iAxsCP].set_ylim(self.CPlim)
axs[iAxsCP].set_ylabel("$c_p$", rotation="horizontal", ha="right", va="center")
# Make legend for viscous vs. inviscid
customLines = [
mpl_lines.Line2D([0], [0], linestyle="-", color="k"),
mpl_lines.Line2D([0], [0], linestyle="--", color="k", linewidth=1.0),
]
axs[iAxsCP].legend(customLines, ["Viscous", "Inviscid"])
# Plot the skin friction coefficient
axs[iAxsCF].plot([min(x) - 1, max(x) + 1], [0.0, 0.0], zorder=-2, color="k", linewidth=0.7, alpha=0.3)
axs[iAxsCF].plot(xCFUpper, CFUpper, color="k", zorder=-1, alpha=0.15)
axs[iAxsCF].plot(xCFLower, CFLower, color="k", zorder=-1, alpha=0.15)
axs[iAxsCF].plot(xCFUpper, CFUpper, color=cpUpColor)
axs[iAxsCF].plot(xCFLower, CFLower, color=cpLowColor)
axs[iAxsCF].set_ylim(self.CFlim)
axs[iAxsCF].set_ylabel("$c_f$", rotation="horizontal", ha="right", va="center")
# Plot the airfoil on the lower axis
axs[iAxsFoil].plot(x, y, color="k", zorder=-1, alpha=0.15)
axs[iAxsFoil].plot(x, y, color=color)
axs[iAxsFoil].set_xlim(self.xlimFoil)
axs[iAxsFoil].set_ylim(self.ylimFoil)
axs[iAxsFoil].set_xlabel("x/c")
axs[iAxsFoil].set_ylabel("y/c", rotation="horizontal", ha="right", va="center")
axs[iAxsFoil].set_aspect("equal")
axs[iAxsFoil].spines["right"].set_visible(False)
axs[iAxsFoil].spines["top"].set_visible(False)
axs[iAxsFoil].yaxis.set_ticks_position("left")
axs[iAxsFoil].xaxis.set_ticks_position("bottom")
if fileName is None and showPlot:
plt.pause(0.5)
self.airfoilFig = fig
self.airfoilAxs = axs
else:
self._updateAirfoilPlot(pause=showPlot)
if fileName is not None:
self.airfoilFig.savefig(fileName)
return self.airfoilFig, self.airfoilAxs
def _updateAirfoilPlot(self, pause=True):
"""
Updates the airfoil plot with current airfoil shape.
Assumes that the current figure is the one with the
airfoil on it and it was the most recently plotted line.
This function should not be called directly by the user
unless you know what you're doing.
Parameters
----------
pause : bool
If true, pauses after updating plot for 0.5 sec. This is necessary
when updating a live plot, but breaks the animation in postprocess.
"""
if plt is None:
raise ImportError("matplotlib is required to use the plotting tools")
# Get data to plot
x = self.coords[:, 0]
y = self.coords[:, 1]
CPUpper = self.sliceData[self.curAP.name]["cp_visc_upper"]
CPUpperInvisc = self.sliceData[self.curAP.name]["cp_invisc_upper"]
xUpper = self.sliceData[self.curAP.name]["x_upper"]
CPLower = self.sliceData[self.curAP.name]["cp_visc_lower"]
CPLowerInvisc = self.sliceData[self.curAP.name]["cp_invisc_lower"]
xLower = self.sliceData[self.curAP.name]["x_lower"]
CFUpper = self.sliceData[self.curAP.name]["cf_upper"]
xCFUpper = self.sliceData[self.curAP.name]["x_cf_upper"]
CFLower = self.sliceData[self.curAP.name]["cf_lower"]
xCFLower = self.sliceData[self.curAP.name]["x_cf_lower"]
iAxsCP = 0
iAxsCF = 1
iAxsFoil = 2
# Compute new plot limits
stackedCP = np.hstack((CPUpper, CPUpperInvisc, CPLower, CPLowerInvisc))
stackedCP = stackedCP[np.isfinite(stackedCP)]
stackedCF = np.hstack((CFUpper, CFLower))
stackedCF = stackedCF[np.isfinite(stackedCF)]
CPlim = [max(stackedCP) + 0.2, min(stackedCP) - 0.2]
CPlim = [max(CPlim[0], self.CPlim[0]), min(CPlim[1], self.CPlim[1])] # inverted axis
CFlim = [min(stackedCF) - 0.002, max(stackedCF) + 0.002]
CFlim = [min(CFlim[0], self.CFlim[0]), max(CFlim[1], self.CFlim[1])]
xlimFoil = [min(x) - 0.01, max(x) + 0.01]
ylimFoil = [min(y) - 0.01, max(y) + 0.01]
xlimFoil = [min(xlimFoil[0], self.xlimFoil[0]), max(xlimFoil[1], self.xlimFoil[1])]
ylimFoil = [min(ylimFoil[0], self.ylimFoil[0]), max(ylimFoil[1], self.ylimFoil[1])]
# CP plot
self.airfoilAxs[iAxsCP].lines[-1].remove()
self.airfoilAxs[iAxsCP].lines[-1].remove()
self.airfoilAxs[iAxsCP].lines[-1].remove()
self.airfoilAxs[iAxsCP].lines[-1].remove()
self.airfoilAxs[iAxsCP].plot(xUpper, CPUpper, color=cpUpColor)
self.airfoilAxs[iAxsCP].plot(xLower, CPLower, color=cpLowColor)
self.airfoilAxs[iAxsCP].plot(xUpper, CPUpperInvisc, "--", color=cpUpColor, linewidth=1.0)
self.airfoilAxs[iAxsCP].plot(xLower, CPLowerInvisc, "--", color=cpLowColor, linewidth=1.0)
self.airfoilAxs[iAxsCP].set_ylim(CPlim)
# CF plot
self.airfoilAxs[iAxsCF].lines[-1].remove()
self.airfoilAxs[iAxsCF].lines[-1].remove()
self.airfoilAxs[iAxsCF].plot(xCFUpper, CFUpper, color=cpUpColor)
self.airfoilAxs[iAxsCF].plot(xCFLower, CFLower, color=cpLowColor)
self.airfoilAxs[iAxsCF].set_ylim(CFlim)
self.airfoilAxs[iAxsFoil].lines[-1].remove()
self.airfoilAxs[iAxsFoil].plot(x, y, color=color)
self.airfoilAxs[iAxsFoil].set_xlim(xlimFoil)
self.airfoilAxs[iAxsFoil].set_ylim(ylimFoil)
if pause:
plt.pause(1e-6)
def _setCoordinates(self, x, y, setComplex=False):
"""
Set the x and y coordinates in the compiled XFOIL layer.
Parameters
----------
x : ndarray (# pts,)
x coordinates of airfoil
y : ndarray (# pts,)
y coordinates of airfoil
setComplex : bool, optional
Set coordinates of the complex version, otherwise will set
the coordinates in the real version, by default False
"""
N = 572 # This is VERY Important: The airfoil input must be a
# FIXED length of 572. Simply set the coordinates up
# to NB and leave the remainder as zeros
NB = len(x)
if setComplex:
dtype = complex
xfoil = self.xfoil_cs
else:
dtype = float
xfoil = self.xfoil
x_input = np.zeros(N, dtype=dtype)
y_input = np.zeros(N, dtype=dtype)
x_input[:NB] = np.array(x).copy()
y_input[:NB] = np.array(y).copy()
xfoil.ci04.nb = NB
xfoil.cr14.xb = x_input
xfoil.cr14.yb = y_input
xfoil.xfoil()
@staticmethod
def _getDefaultOptions():
return {
"maxIters": [int, 100], # maximum iterations for XFOIL solver
"writeCoordinates": [bool, True], # whether to write coordinates to .dat file when `writeSolution` called
"writeSliceFile": [bool, True], # whether or not to save chordwise data
"writeSolution": [bool, False], # whether or not to call writeSolution after each call to XFOIL
"plotAirfoil": [bool, False], # whether to plot airfoil while running
"outputDirectory": [str, "."], # where to put output files
"numberSolutions": [bool, True], # whether to add call counter to output file names
"xTrip": [
np.ndarray,
np.full(2, np.NaN),
], # boundary layer trip x coordinate of upper and lower surface, respectively (two-element array)
}
@staticmethod
def _readDat(filename, headerlines=0):
"""
Reads a '.dat' style airfoil coordinate file,
with each coordinate on a new line and each
line containing an xy pair separate by whitespace.
Parameters
----------
filename : str
dat file from which to read data
headerlines : int
the number of lines to skip at the beginning of the file to reach the coordinates
Returns
-------
X : Ndarray [N,2]
The coordinates read from the file
"""
with open(filename, "r") as f:
for _i in range(headerlines):
f.readline()
r = []
while True:
line = f.readline()
if not line:
break # end of file
if line.isspace():
break # blank line
r.append([float(s) for s in line.split()])
X = np.array(r)
return X