diff --git a/baseclasses/problems/pyAero_problem.py b/baseclasses/problems/pyAero_problem.py index 2489049..e93e758 100644 --- a/baseclasses/problems/pyAero_problem.py +++ b/baseclasses/problems/pyAero_problem.py @@ -10,7 +10,7 @@ from .ICAOAtmosphere import ICAOAtmosphere from .FluidProperties import FluidProperties from ..utils import CaseInsensitiveDict, Error - +from collections import defaultdict class AeroProblem(FluidProperties): """ @@ -330,9 +330,32 @@ def __init__(self, name, **kwargs): if getattr(self, var) is not None: self.possibleDVs.add(var) - BCVarFuncs = ["Pressure", "PressureStagnation", "Temperature", "TemperatureStagnation", "Thrust", "Heat"] + BCVarFuncs = [ + "Temperature", + "Pressure", + "Density", + "TemperatureStagnation", + "PressureStagnation", + "DensityStagnation", + "VelocityX", + "VelocityY", + "VelocityZ", + "VelocityR", + "VelocityTheta", + "VelocityUnitVectorX", + "VelocityUnitVectorY", + "VelocityUnitVectorZ", + "VelocityUnitVectorR", + "VelocityUnitVectorTheta", + "VelocityAngleX", + "VelocityAngleY", + "VelocityAngleZ", + ] self.possibleBCDVs = set(BCVarFuncs) + actuatorFuncs = ["Thrust", "Torque", "Heat"] + self.possibleActuatorDVs = set(actuatorFuncs) + # Now determine the possible functions. Any possible design # variable CAN also be a function (pass through) self.possibleFunctions = set(self.possibleDVs) @@ -351,7 +374,8 @@ def __init__(self, name, **kwargs): # Storage of BC varible values # vars are keyed by (bcVarName, Family) - self.bcVarData = {} + self.bc_data = {} + self.actuatorData = {} def _setStates(self, inputDict): """ @@ -459,13 +483,67 @@ def _setStates(self, inputDict): "to correctly specify the aerodynamic state" ) - def setBCVar(self, varName, value, familyName): + def setBCVar(self, varName, value, groupName): """ set the value of a BC variable on a specific variable """ + if not groupName in self.bc_data.keys(): + self.bc_data[groupName] = {} + + self.bc_data[groupName][varName] = value + + def getBCVars(self): + return self.bc_data + + # def setBCArray(self, groupName, varName, dataArray, patch=None): + # if patch == None: + # # assume the data is set on the first patch in this group + # patches = self.BCData[groupName][varName].keys() + # self.BCData[groupName][varName][patches[0]] = dataArray + # else: + # self.BCData[groupName][varName][patch] = dataArray + + def setActuatorVar(self, varName, value, groupName): + if not groupName in self.actuatorData.keys(): + self.actuatorData[groupName] = {} + + self.actuatorData[groupName][varName] = value + + def getActuatorVars(self): + return self.actuatorData + + # def addBCDV(self, bc_var, value=None, lower=None, upper=None, scale=1.0, + # name=None, offset=0.0, dvOffset=0.0, addToPyOpt=True, familyGroup=None, + # units=None, **kwargs): + # """ + # add a set of bcdata was a design variable + # """ + + # if bc_var not in self.possibleBCDVs: + # raise ValueError('%s is not a valid design variable' % key) + + # # get a numpy array by appliying thee filters to the bc data + # data = self.bc_data.getBCArraysFlatData( bc_var, familyGroup=familyGroup) + + # if name is None: + # dvName = '%s_%s_%s' % (key, familyGroup, self.name) + # else: + # dvName = name - self.bcVarData[varName, familyName] = value - print("update bc", value) + # # combine the data + + # if value is None: + # value = data + # else: + # # insure that the given size is correct if it is not a float + # # + # if (not isinstance(value, float)) and value.size == data.size: + # raise Error('give value for {} is not ht right size, {} given {} expected'.format(bc_var, value.size, data.size)) + + # dv = aeroDV(bc_var, value, lower, upper, scale, offset, + # dvOffset, addToPyOpt, familyGroup, units) + # self.DVs[dvName] = dv + # return def addDV( self, @@ -478,7 +556,7 @@ def addDV( offset=0.0, dvOffset=0.0, addToPyOpt=True, - family=None, + familyGroup=None, units=None, ): """ @@ -548,7 +626,7 @@ def addDV( >>> ap.addDV('alpha', value=2.5, lower=0.0, upper=10.0, scale=0.1) """ - if (key not in self.allVarFuncs) and (key not in self.possibleBCDVs): + if (key not in self.allVarFuncs) and (key not in self.possibleBCDVs) and (key not in self.possibleActuatorDVs): raise ValueError("%s is not a valid design variable" % key) # First check if we are allowed to add the DV: @@ -560,19 +638,71 @@ def addDV( "(...,alpha=value, ...) must be given. The list of " "possible DVs are: %s." % (key, repr(self.possibleDVs)) ) - if key in self.possibleBCDVs: - if family is None: - raise Error("The family must be given for BC design variables") + + if key in self.possibleBCDVs or key in self.possibleActuatorDVs: + if familyGroup is None: + raise Error( + "The familyGroup must be given for BC or actuator\ + design variables" + ) if name is None: - dvName = f"{key}_{family}_{self.name}" + dvName = f"{key}_{familyGroup}_{self.name}" else: dvName = name if value is None: - if (key, family) not in self.bcVarData: - raise Error("The value must be given or set using the setBCVar routine") - value = self.bcVarData[key, family] + if key in self.possibleBCDVs: + try: + value = self.bc_data[familyGroup][key] + except KeyError: + raise Error("The value must be given or set using the setBCVar routine") + else: + try: + value = self.actuatorData[familyGroup][key] + except KeyError: + raise Error("The value must be given or set using the setActuatorVar routine") + + # # the value of the BCData[familyGroup][key] maybe a dictionary of data for each patch + # # to accomadate this we will add a variable for each entry in the dictionary + # if isinstance(value, dict) and key in self.possibleBCDVs: + + # # use a little bit o recursion + # for dict_key in value: + # dict_key_name = str(dict_key).replace(" ", "") + # dict_dvName = '%s_%s' % (dvName , dict_key_name) + + # if isinstance(lower, dict): + # lower_val = lower[dict_key] + # else: + # lower_val = lower + # if isinstance(upper, dict): + # upper_val = upper[dict_key] + # else: + # upper_val = upper + # if isinstance(scale, dict): + # scale_val = scale[dict_key] + # else: + # scale_val = scale + # if isinstance(offset, dict): + # offset_val = offset[dict_key] + # else: + # offset_val = offset + # if isinstance(dvOffset, dict): + # dvOffset_val = dvOffset[dict_key] + # else: + # dvOffset_val = dvOffset + + # # self.addDV(key, , lower_val, upper_val, scale_val, dvName, + # # offset_val, dvOffset_val, addToPyOpt, familyGroup, + # # units) + + # self.DVs[dict_dvName] = aeroDV(key, value[dict_key], lower_val, upper_val, scale_val, offset_val, + # dvOffset_val, addToPyOpt, familyGroup, units) + # self.DVs[dict_dvName].dict_key = dict_key + + # return + else: if name is None: dvName = key + "_%s" % self.name @@ -581,9 +711,9 @@ def addDV( if value is None: value = getattr(self, key) - family = None + familyGroup = None - self.DVs[dvName] = aeroDV(key, value, lower, upper, scale, offset, dvOffset, addToPyOpt, family, units) + self.DVs[dvName] = aeroDV(key, value, lower, upper, scale, offset, dvOffset, addToPyOpt, familyGroup, units) def updateInternalDVs(self): """ @@ -609,18 +739,24 @@ def setDesignVars(self, x): for dvName in self.DVs: if dvName in x: - key = self.DVs[dvName].key - family = self.DVs[dvName].family - value = x[dvName] + self.DVs[dvName].offset - if family is None: - setattr(self, key, value) + dv = self.DVs[dvName] + key = dv.key + family = dv.family + value = x[dvName] + dv.offset + if key in self.possibleBCDVs: + self.bc_data[family][key] = value + + elif key in self.possibleActuatorDVs: + self.actuatorData[family][key] = value else: - self.bcVarData[key, family] = value + setattr(self, key, value) - try: # To set in the DV as well if the DV exists: - self.DVs[dvName].value = x[dvName] - except: # noqa - pass # DV doesn't exist + dv.value = x[dvName] + + # try: # To set in the DV as well if the DV exists: + # self.DVs[dvName].value = x[dvName] + # except: # noqa + # pass # DV doesn't exist def addVariablesPyOpt(self, optProb): """ @@ -824,50 +960,6 @@ def altitude(self, value): self._setStates({"altitude": value}) self._set_aeroDV_val("altitude", value) - # def _update(self): - # """ - # Try to finish the complete state: - # """ - - # if self.T is not None: - # self.a = numpy.sqrt(self.gamma*self.R*self.T) - # if self.englishUnits: - # mu = (self.muSuthDim * ( - # (self.TSuthDim + self.SSuthDim) / (self.T/1.8 + self.SSuthDim)) * - # (((self.T/1.8)/self.TSuthDim)**1.5)) - # self.mu = mu / 47.9 - # else: - # self.mu = (self.muSuthDim * ( - # (self.TSuthDim + self.SSuthDim) / (self.T + self.SSuthDim)) * - # ((self.T/self.TSuthDim)**1.5)) - - # if self.mach is not None and self.a is not None: - # self.V = self.mach * self.a - - # if self.a is not None and self.V is not None: - # self.__dict__['mach'] = self.V/self.a - - # if self.P is not None and self.T is not None: - # self.__dict__['rho'] = self.P/(self.R*self.T) - - # if self.rho is not None and self.T is not None: - # self.__dict__['P'] = self.rho*self.R*self.T - - # if self.rho is not None and self.P is not None: - # self.__dict__['T'] = self.P /(self.rho*self.R) - - # if self.mu is not None and self.rho is not None: - # self.nu = self.mu / self.rho - - # if self.rho is not None and self.V is not None: - # self.q = 0.5*self.rho*self.V**2 - - # if self.rho is not None and self.V is not None and self.mu is not None: - # self.__dict__['re'] = self.rho*self.V/self.mu - - # if self.re is not None and self.mu is not None and self.V is not None: - # self.__dict__['rho'] = self.re*self.mu/self.V - def _updateFromRe(self): """ update the full set of states from M,T,P @@ -957,6 +1049,96 @@ def _getDVSens(self, func): return rDict + def evalDVsSensFwd(self, xDvDot): + """ + fwd mod deriv of setDesignVars + """ + aeroDvsDot = {} + bcDvsDot = defaultdict(lambda: {}) + actDvsDot = defaultdict(lambda: {}) + + for dvName in self.DVs: + if dvName in xDvDot: + dv = self.DVs[dvName] + key = dv.key + family = dv.family + value = xDvDot[dvName] + + if key in self.possibleBCDVs: + bcDvsDot[family][key] = value + + elif key in self.possibleActuatorDVs: + actDvsDot[family][key] = value + else: + aeroDvsDot[key] = value + + return aeroDvsDot, bcDvsDot, actDvsDot + + def evalDVsSensBwd(self, aeroDvBar, BCDataBar, actDataBar): + """This internal furncion is used to convert the raw array output from + the matrix-free product bwd routine into the required + dictionary format.""" + + DVbar = {} + for dvName in self.DVs: + key = self.DVs[dvName].key.lower() + + tmp = {} + if key == "altitude": + # This design variable is special. It combines changes + # in temperature, pressure and density into a single + # variable. Since we have derivatives for T, P and + # rho, we simply chain rule it back to the the + # altitude variable. + self.evalFunctionsSens(tmp, ["P", "T", "rho"]) + + # Extract the derivatives wrt the independent + # parameters in ADflow + dIdP = aeroDvBar["p"] + dIdT = aeroDvBar["t"] + dIdrho = aeroDvBar["rho"] + + # Chain-rule to get the final derivative: + DVbar[dvName] = ( + tmp[self["P"]][dvName] * dIdP + tmp[self["T"]][dvName] * dIdT + tmp[self["rho"]][dvName] * dIdrho + ) + + elif key == "mach": + self.evalFunctionsSens(tmp, ["P", "rho"]) + # Simular story for Mach: It is technically possible + # to use a mach number for a fixed RE simulation. For + # the RE to stay fixed and change the mach number, the + # 'P' and 'rho' must also change. We have to chain run + # this dependence back through to the final mach + # derivative. When Mach number is used with altitude + # or P and T, this calc is unnecessary, but won't do + # any harm. + dIdP = aeroDvBar["p"] + dIdrho = aeroDvBar["rho"] + + # Chain-rule to get the final derivative: + DVbar[dvName] = tmp[self["P"]][dvName] * dIdP + tmp[self["rho"]][dvName] * dIdrho + aeroDvBar["mach"] + # if the variable is an BC DV, get the data from the BCDataBar + # dict for each patch + elif key in [k.lower() for k in self.possibleBCDVs]: + fam = self.DVs[dvName].family + dv_key = self.DVs[dvName].key + + DVbar[dvName] = BCDataBar[fam][dv_key] + + # if the variable is an actuator DV, get the data from the actDataBar + elif key in [k.lower() for k in self.possibleActuatorDVs]: + fam = self.DVs[dvName].family + dv_key = self.DVs[dvName].key + DVbar[dvName] = actDataBar[fam][dv_key] + + else: + DVbar[dvName] = aeroDvBar[key] + + return DVbar + + # return DVbar + class aeroDV: """ @@ -974,3 +1156,18 @@ def __init__(self, key, value, lower, upper, scale, offset, dvOffset, addToPyOpt self.addToPyOpt = addToPyOpt self.family = family self.units = units + + if isinstance(value, float): + self.size = 1 + elif isinstance(value, numpy.ndarray): + self.size = value.size + else: + size = None + # TODO raise an error? + + def __repr__(self): + """Returns representation of the object""" + msg = "{}({!r}) {}".format(self.__class__.__name__, self.key, self.value) + if self.units is not None: + msg += self.units + return msg \ No newline at end of file diff --git a/baseclasses/solvers/pyAero_solver.py b/baseclasses/solvers/pyAero_solver.py index 26ad7a4..731083b 100644 --- a/baseclasses/solvers/pyAero_solver.py +++ b/baseclasses/solvers/pyAero_solver.py @@ -41,6 +41,10 @@ def __init__( """ AeroSolver Class Initialization """ + + self.families = CaseInsensitiveDict() + self.familyGroups = CaseInsensitiveDict() + # Setup option info super().__init__( name, @@ -284,10 +288,25 @@ def resetFlow(self): pass - def addFamilyGroup(self, groupName, families): + def addFamilyGroup(self, groupName, groups): """Add a custom grouping of families called groupName. The groupName - must be distinct from the existing families. All families must - in the 'families' list must be present in the CGNS file. + must be distinct from the existing family groups names. + # All group in the 'families' list must be present in the CGNS file. + + + + ---------------------------------- + self.families | 'inlet' 'outlet' 'nozzle_wall' | + (from surface definition) | [1] [2] [3] | + ---------------------------------- + + ----------------------------------------------------------------- + self.familyGroups | 'inlet' 'outlet' 'nozzle_wall' 'allwalls' 'allsurfaces' | + (collection of | [1] [2] [3] [3] [1, 2, 3] | + families) ----------------------------------------------------------------- + + all familes are family groups since each is a collection of a single + family (itself) Parameters ---------- @@ -306,23 +325,23 @@ def addFamilyGroup(self, groupName, families): # We can actually allow for nested groups. That is, an entry # in families may already be a group added in a previous call. indices = [] - for fam in families: - if fam not in self.families: - raise Error( - "The specified family '%s' for group '%s', does " - "not exist in the cgns file or has " - "not already been added. The current list of " - "families (original and grouped) is: %s" % (fam, groupName, repr(self.families.keys())) - ) + for group in groups: + if group.lower() not in self.familyGroups: + raise Error("The specified familyGroup '%s' for group '%s', does " + "not exist in the cgns file or has " + "not already been added. The current list of " + "families is: %s and the current list of family" + "groups is: %s"%(group, groupName, repr(self.families.keys()), + repr(self.familyGroups.keys())) ) - indices.extend(self.families[fam]) + indices.extend(self.familyGroups[group]) # It is very important that the list of families is sorted # becuase in fortran we always use a binary search to check if # a famID is in the list. - self.families[groupName] = sorted(numpy.unique(indices)) + self.familyGroups[groupName] = sorted(numpy.unique(indices)) - def getSurfaceCoordinates(self, group_name): + def getSurfaceCoordinates(self,group_name): """ Return the set of surface coordinates cooresponding to a Particular group name @@ -473,6 +492,18 @@ def printFamilyList(self): pprint(self.families) + def printFamilyGroupList(self): + """ + Print a nicely formatted dictionary of the family names + """ + from pprint import pprint + pprint(self.familyGroups) + + +# -------------------------- +# Private Utility functions +# -------------------------- + # -------------------------- # Private Utility functions # -------------------------- @@ -482,10 +513,19 @@ def _getFamilyList(self, groupName): if groupName is None: groupName = self.allFamilies - if groupName not in self.families: - raise Error( - "'%s' is not a family in the CGNS file or has not been added" - " as a combination of families" % groupName - ) + if groupName not in self.familyGroups: + raise Error("'%s' is not a family in the CGNS file or has not been added" + " as a combination of families"%groupName) + + return self.familyGroups[groupName] + + def _getFamiliesInFamilyGroup(self, groupName): + famlist = self._getFamilyList(groupName) + + groupFamilies = [] + for famID in famlist: + for family in self.families: + if famID == self.families[family]: + groupFamilies.append(family) - return self.families[groupName] + return groupFamilies diff --git a/baseclasses/testing/pyRegTest.py b/baseclasses/testing/pyRegTest.py index 9efcce7..659e25e 100644 --- a/baseclasses/testing/pyRegTest.py +++ b/baseclasses/testing/pyRegTest.py @@ -332,16 +332,17 @@ def _add_dict(self, dict_name, d, full_name, db=None, **kwargs): db[dict_name] = {} elif dict_name not in db.keys(): raise ValueError(f"The key '{dict_name}' was not found in the reference file!") - + for key in sorted(d.keys()): - full_name = f"{full_name}: {key}" + full_name_key = f"{full_name}: {key}" + if isinstance(d[key], bool): - self._add_values(key, int(d[key]), rtol=rtol, atol=atol, db=db[dict_name], full_name=full_name) + self._add_values(key, int(d[key]), rtol=rtol, atol=atol, db=db[dict_name], full_name=full_name_key) elif isinstance(d[key], dict): # do some good ol' fashion recursion - self._add_dict(key, d[key], full_name, rtol=rtol, atol=atol, db=db[dict_name]) + self._add_dict(key, d[key], full_name_key, rtol=rtol, atol=atol, db=db[dict_name]) else: - self._add_values(key, d[key], rtol=rtol, atol=atol, db=db[dict_name], full_name=full_name) + self._add_values(key, d[key], rtol=rtol, atol=atol, db=db[dict_name], full_name=full_name_key) # This strategy of dealing with error propagation to multiple procs is taken directly form openMDAO.utils; diff --git a/baseclasses/utils/__init__.py b/baseclasses/utils/__init__.py index b2bd75e..fd81351 100644 --- a/baseclasses/utils/__init__.py +++ b/baseclasses/utils/__init__.py @@ -1,6 +1,6 @@ from .containers import CaseInsensitiveSet, CaseInsensitiveDict from .error import Error -from .utils import getPy3SafeString, pp +from .utils import getPy3SafeString, pp, printArgs from .fileIO import writeJSON, readJSON, writePickle, readPickle, redirectIO, redirectingIO __all__ = [ @@ -15,4 +15,5 @@ "readPickle", "redirectIO", "redirectingIO", + "printArgs", ] diff --git a/baseclasses/utils/utils.py b/baseclasses/utils/utils.py index 2506521..f1c9df9 100644 --- a/baseclasses/utils/utils.py +++ b/baseclasses/utils/utils.py @@ -41,3 +41,31 @@ def pp(obj, comm=None, flush=True): # we use pformat to get the string and then call print manually, that way we can flush if we need to pprint_str = pformat(obj) print(pprint_str, flush=flush) + +def printArgs(args): + """ + Prints the arguments passed to the script. + + Parameters + ---------- + args : argparse.Namespace + The object holding the arguments passed to the script from + args = parser.parse_args() + """ + args_dict = vars(args) + longest_key = max(args_dict.keys(), key=len) + longest_val = max([str(x) for x in args_dict.values()], key=len) + divider = " : " + + box_width = len(longest_key) + len(longest_val) + len(divider) + + bar = "-" * box_width + + # add title in the middle of the bar to create header + title = " Arguments " + header = bar[:(len(bar) - len(title))//2] + title + bar[(len(bar)+len(title))//2:] + + print(header) + for arg, arg_val in args_dict.items(): + print(f"{arg:{len(longest_key)}}" + divider + f"{arg_val}") + print(bar)