Skip to content

Commit

Permalink
Fixed optFuncs.
Browse files Browse the repository at this point in the history
  • Loading branch information
@friedenhe
friedenhe committed Jan 13, 2024
1 parent f713943 commit bd006d9
Showing 1 changed file with 59 additions and 32 deletions.
91 changes: 59 additions & 32 deletions dafoam/optFuncs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,11 +42,13 @@ def calcObjFuncValues(xDV):
funcs = {}

# Set the current design variables in the DV object
DVGeo.setDesignVars(xDV)
if DVGeo is not None:
DVGeo.setDesignVars(xDV)
DASolver.setInternalDesignVars(xDV)

# Evaluate the geometric constraints and add them to the funcs dictionary
DVCon.evalFunctions(funcs)
if DVCon is not None:
DVCon.evalFunctions(funcs)

# Solve the CFD problem
DASolver()
Expand Down Expand Up @@ -87,7 +89,8 @@ def calcObjFuncValuesMP(xDV):
funcsMP = {}

# Set the current design variables in the DV object
DVGeo.setDesignVars(xDV)
if DVGeo is not None:
DVGeo.setDesignVars(xDV)
DASolver.setInternalDesignVars(xDV)

nMultiPoints = DASolver.getOption("nMultiPoints")
Expand All @@ -99,7 +102,8 @@ def calcObjFuncValuesMP(xDV):
funcs = {}

# Evaluate the geometric constraints and add them to the funcs dictionary
DVCon.evalFunctions(funcs)
if DVCon is not None:
DVCon.evalFunctions(funcs)

# set the multi point condition provided by users in the
# runScript.py script. This function should define what
Expand Down Expand Up @@ -156,11 +160,13 @@ def calcObjFuncValuesUnsteady(xDV):
funcs = {}

# Set the current design variables in the DV object
DVGeo.setDesignVars(xDV)
if DVGeo is not None:
DVGeo.setDesignVars(xDV)
DASolver.setInternalDesignVars(xDV)

# Evaluate the geometric constraints and add them to the funcs dictionary
DVCon.evalFunctions(funcs)
if DVCon is not None:
DVCon.evalFunctions(funcs)

# Solve the CFD problem
DASolver()
Expand Down Expand Up @@ -206,7 +212,8 @@ def calcObjFuncSens(xDV, funcs):
funcsSens = {}

# Evaluate the geometric constraint derivatives
DVCon.evalFunctionsSens(funcsSens)
if DVCon is not None:
DVCon.evalFunctionsSens(funcsSens)

# Solve the adjoint
if DASolver.getOption("unsteadyAdjoint")["mode"] == "timeAccurate":
Expand Down Expand Up @@ -266,7 +273,8 @@ def calcObjFuncSensMP(xDV, funcs):
funcsSens = {}

# Evaluate the geometric constraint derivatives
DVCon.evalFunctionsSens(funcsSens)
if DVCon is not None:
DVCon.evalFunctionsSens(funcsSens)

# set the state vector for case i
DASolver.setMultiPointField(i)
Expand Down Expand Up @@ -356,7 +364,8 @@ def calcObjFuncSensUnsteady(xDV, funcs):
funcsSens = {}

# Evaluate the geometric constraint derivatives
DVCon.evalFunctionsSens(funcsSens)
if DVCon is not None:
DVCon.evalFunctionsSens(funcsSens)

# set the state vector for case i
DASolver.setTimeInstanceField(i)
Expand Down Expand Up @@ -397,9 +406,14 @@ def runPrimal(objFun=calcObjFuncValues):
Just run the primal
"""

xDV = DVGeo.getValues()
xDV = {}
if DVGeo is not None:
xDV = DVGeo.getValues()
iDV = DASolver.getInternalDVDict()
allDV = {**xDV, **iDV}

funcs = {}
funcs, fail = objFun(xDV)
funcs, fail = objFun(allDV)

return funcs, fail

Expand All @@ -410,18 +424,23 @@ def runAdjoint(objFun=calcObjFuncValues, sensFun=calcObjFuncSens, fileName=None)
"""

DASolver.runColoring()
xDV = DVGeo.getValues()
xDV = {}
if DVGeo is not None:
xDV = DVGeo.getValues()
iDV = DASolver.getInternalDVDict()
allDV = {**xDV, **iDV}

funcs = {}
funcs, fail = objFun(xDV)
funcs, fail = objFun(allDV)
funcsSens = {}
funcsSens, fail = sensFun(xDV, funcs)
funcsSens, fail = sensFun(allDV, funcs)

# Optionally, we can write the sensitivity to a file if fileName is provided
if fileName is not None:
if gcomm.rank == 0:
fOut = open(fileName, "w")
for funcName in evalFuncs:
for shapeVar in xDV:
for shapeVar in allDV:
fOut.write(funcName + " " + shapeVar + "\n")
try:
nDVs = len(funcsSens[funcName][shapeVar])
Expand Down Expand Up @@ -463,24 +482,28 @@ def solveCL(CL_star, alphaName, liftName, objFun=calcObjFuncValues, eps=1e-2, to
Info("+--------------------------------------------------------------------------+")
Info("eps: %g tol: %g maxit: %g" % (eps, tol, maxit))

xDVs = DVGeo.getValues()
alpha = xDVs[alphaName]
xDV = {}
if DVGeo is not None:
xDV = DVGeo.getValues()
iDV = DASolver.getInternalDVDict()
allDV = {**xDV, **iDV}
alpha = allDV[alphaName]

for i in range(maxit):
# Solve the CFD problem
xDVs[alphaName] = alpha
allDV[alphaName] = alpha
funcs = {}
funcs, fail = objFun(xDVs)
funcs, fail = objFun(allDV)
CL0 = funcs[liftName]
Info("alpha: %f, CL: %f" % (alpha.real, CL0))
if abs(CL0 - CL_star) / CL_star < tol:
Info("Completed! alpha = %f" % alpha.real)
return alpha.real
# compute sens
alphaVal = alpha + eps
xDVs[alphaName] = alphaVal
allDV[alphaName] = alphaVal
funcsP = {}
funcsP, fail = objFun(xDVs)
funcsP, fail = objFun(allDV)
CLP = funcsP[liftName]
deltaAlpha = (CL_star - CL0) * eps / (CLP - CL0)
alpha += deltaAlpha
Expand All @@ -493,36 +516,40 @@ def calcFDSens(objFun=calcObjFuncValues, fileName=None):
Compute finite-difference sensitivity
"""

xDV = DVGeo.getValues()
xDV = {}
if DVGeo is not None:
xDV = DVGeo.getValues()
iDV = DASolver.getInternalDVDict()
allDV = {**xDV, **iDV}

# gradFD
deltaX = DASolver.getOption("adjPartDerivFDStep")["FFD"]
# initialize gradFD
gradFD = {}
for funcName in evalFuncs:
gradFD[funcName] = {}
for shapeVar in xDV:
for shapeVar in allDV:
try:
nDVs = len(xDV[shapeVar])
nDVs = len(allDV[shapeVar])
except Exception:
nDVs = 1
gradFD[funcName][shapeVar] = np.zeros(nDVs)
if gcomm.rank == 0:
print("-------FD----------", deltaX, flush=True)

for shapeVar in xDV:
for shapeVar in allDV:
try:
nDVs = len(xDV[shapeVar])
nDVs = len(allDV[shapeVar])
except Exception:
nDVs = 1
for i in range(nDVs):
funcp = {}
funcm = {}
xDV[shapeVar][i] += deltaX
funcp, fail = objFun(xDV)
xDV[shapeVar][i] -= 2.0 * deltaX
funcm, fail = objFun(xDV)
xDV[shapeVar][i] += deltaX
allDV[shapeVar][i] += deltaX
funcp, fail = objFun(allDV)
allDV[shapeVar][i] -= 2.0 * deltaX
funcm, fail = objFun(allDV)
allDV[shapeVar][i] += deltaX
for funcName in evalFuncs:
gradFD[funcName][shapeVar][i] = (funcp[funcName] - funcm[funcName]) / (2.0 * deltaX)
Info(gradFD)
Expand All @@ -532,7 +559,7 @@ def calcFDSens(objFun=calcObjFuncValues, fileName=None):
fOut = open(fileName, "w")
fOut.write("DeltaX: " + str(deltaX) + "\n")
for funcName in evalFuncs:
for shapeVar in xDV:
for shapeVar in allDV:
fOut.write(funcName + " " + shapeVar + "\n")
try:
nDVs = len(gradFD[funcName][shapeVar])
Expand Down

0 comments on commit bd006d9

Please sign in to comment.