Added internalDV (#561)

* Added default val for DARegModel.

* Added internal design vars.

* Fixed an issue.

* Fixed an issue and add new tests for internalDV.

* Fixed some issues for internalDV.

dafoam/optFuncs.py

@@ -43,7 +43,7 @@ def calcObjFuncValues(xDV):
 
     # Set the current design variables in the DV object
     DVGeo.setDesignVars(xDV)
-    DASolver.setDesignVars(xDV)
+    DASolver.setInternalDesignVars(xDV)
 
     # Evaluate the geometric constraints and add them to the funcs dictionary
     DVCon.evalFunctions(funcs)
@@ -88,7 +88,7 @@ def calcObjFuncValuesMP(xDV):
 
     # Set the current design variables in the DV object
     DVGeo.setDesignVars(xDV)
-    DASolver.setDesignVars(xDV)
+    DASolver.setInternalDesignVars(xDV)
 
     nMultiPoints = DASolver.getOption("nMultiPoints")
 
@@ -157,7 +157,7 @@ def calcObjFuncValuesUnsteady(xDV):
 
     # Set the current design variables in the DV object
     DVGeo.setDesignVars(xDV)
-    DASolver.setDesignVars(xDV)
+    DASolver.setInternalDesignVars(xDV)
 
     # Evaluate the geometric constraints and add them to the funcs dictionary
     DVCon.evalFunctions(funcs)

dafoam/pyDAFoam.py

@@ -523,6 +523,7 @@ def __init__(self):
             "outputLowerBound": -1e8,
             "activationFunction": "sigmoid",
             "printInputRange": False,
+            "defaultOutputValue": 0.0,
         }
 
         # *********************************************************************************************
@@ -887,6 +888,9 @@ def __init__(self, comm=None, options=None):
         # initialize the dRdWTPC dict for unsteady adjoint
         self.dRdWTPCUnsteady = None
 
+        # initialize the user defined internal dvs dict
+        self.internalDV = {}
+
         if self.getOption("tensorflow")["active"]:
             TensorFlowHelper.options = self.getOption("tensorflow")
             TensorFlowHelper.initialize()
@@ -958,6 +962,9 @@ def __call__(self):
                 if dvType == "FFD":
                     self.calcFFD2XvSeedVec()
 
+        # now call the internal design var function to update DASolver parameters
+        self.runInternalDVFunc()
+
         # update the primal boundary condition right before calling solvePrimal
         self.setPrimalBoundaryConditions()
 
@@ -1486,15 +1493,23 @@ def evalFunctionsSens(self, funcsSens, evalFuncs=None):
         >>> CFDsolver.evalFunctionsSens(funcSens, ['CD', 'CL'])
         """
 
-        if self.DVGeo is None:
-            raise Error("DVGeo not set!")
+        for funcName in evalFuncs:
+            funcsSens[funcName] = {}
 
-        dvs = self.DVGeo.getValues()
+        if self.DVGeo is not None:
+
+            dvs = self.DVGeo.getValues()
+
+            for funcName in evalFuncs:
+                for dvName in dvs:
+                    nDVs = len(dvs[dvName])
+                    funcsSens[funcName][dvName] = np.zeros(nDVs, self.dtype)
+                    for i in range(nDVs):
+                        funcsSens[funcName][dvName][i] = self.adjTotalDeriv[funcName][dvName][i]
 
         for funcName in evalFuncs:
-            funcsSens[funcName] = {}
-            for dvName in dvs:
-                nDVs = len(dvs[dvName])
+            for dvName in self.internalDV:
+                nDVs = len(self.internalDV[dvName]["init"])
                 funcsSens[funcName][dvName] = np.zeros(nDVs, self.dtype)
                 for i in range(nDVs):
                     funcsSens[funcName][dvName][i] = self.adjTotalDeriv[funcName][dvName][i]
@@ -1525,6 +1540,72 @@ def setDVGeo(self, DVGeo):
 
         self.DVGeo = DVGeo
 
+    def setInternalDesignVars(self, xDVs):
+        """
+        Set the internal design variables.
+        """
+
+        for dvName in self.internalDV:
+            for i in range(len(self.internalDV[dvName]["value"])):
+                self.internalDV[dvName]["value"][i] = xDVs[dvName][i]
+
+    def getInternalDVDict(self):
+        """
+        Get the internal design variable values
+        """
+        internalDVDict = {}
+        for dvName in self.internalDV:
+            internalDVDict[dvName] = self.internalDV[dvName]["value"]
+
+        return internalDVDict
+
+    def addInternalDV(self, dvName, dvInit, dvFunc, lower, upper, scale):
+        """
+        Add design variable.
+
+        Inputs:
+            dvName: [str] Name of of the design variable
+            dvInit: [array] Initial values for the design variables
+            fvFunc: [function] A function that define how to apply the change
+                     based on the design variable values. The form of this func
+                     is dvFunc(dvVal, DASolver)
+            lower/upper [scalar] The lower/upper bound of the DV
+            scale [scalar] The scaling factor for the DV
+        """
+        self.internalDV[dvName] = {}
+        self.internalDV[dvName]["init"] = dvInit
+        self.internalDV[dvName]["func"] = dvFunc
+        self.internalDV[dvName]["lower"] = lower
+        self.internalDV[dvName]["upper"] = upper
+        self.internalDV[dvName]["scale"] = scale
+        nInternalDVs = len(dvInit)
+        self.internalDV[dvName]["value"] = np.zeros(nInternalDVs)
+        for i in range(nInternalDVs):
+            self.internalDV[dvName]["value"][i] = self.internalDV[dvName]["init"][i]
+
+    def runInternalDVFunc(self):
+        """
+        call the design variable function
+        """
+        for dvName in self.internalDV:
+            Info("Calling internal design variable functioins %s" % dvName)
+            func = self.internalDV[dvName]["func"]
+            dvVal = self.internalDV[dvName]["value"]
+            func(dvVal, self)
+
+    def addVariablesPyOpt(self, optProb):
+        """
+        Add the design variable for optProb. This is similar to the
+        function in DVGeo
+        """
+        for dvName in self.internalDV:
+            dvInit = self.internalDV[dvName]["init"]
+            lower = self.internalDV[dvName]["lower"]
+            upper = self.internalDV[dvName]["upper"]
+            scale = self.internalDV[dvName]["scale"]
+            nDVs = len(dvInit)
+            optProb.addVarGroup(dvName, nDVs, "c", value=dvInit, lower=lower, upper=upper, scale=scale)
+
     def addFamilyGroup(self, groupName, families):
         """
         Add a custom grouping of families called groupName. The groupName
@@ -1714,15 +1795,6 @@ def printFamilyList(self):
         """
         Info(self.families)
 
-    def setDesignVars(self, x):
-        """
-        Set the internal design variables.
-        At the moment we don't have any internal DVs to set.
-        """
-        pass
-
-        return
-
     def _initializeOptions(self, options):
         """
         Initialize the options passed into pyDAFoam
@@ -2174,8 +2246,16 @@ def calcTotalDerivsFFD(self, objFuncName, designVarName, dFScaling=1.0, accumula
 
     def calcTotalDerivsField(self, objFuncName, designVarName, fieldType, dFScaling=1.0, accumulateTotal=False):
 
-        xDV = self.DVGeo.getValues()
-        nDVs = len(xDV[designVarName])
+        nDVs = 0
+        if self.DVGeo is not None:
+            xDV = self.DVGeo.getValues()
+            if designVarName in xDV:
+                nDVs = len(xDV[designVarName])
+        elif designVarName in self.internalDV:
+            nDVs = len(self.internalDV[designVarName]["init"])
+        else:
+            raise Error("design variable %s not found..." % designVarName)
+
         if fieldType == "scalar":
             fieldComp = 1
         elif fieldType == "vector":
@@ -2290,8 +2370,17 @@ def calcTotalDerivsACT(self, objFuncName, designVarName, designVarType, dFScalin
 
     def calcTotalDerivsRegPar(self, objFuncName, designVarName, dFScaling=1.0, accumulateTotal=False):
 
-        xDV = self.DVGeo.getValues()
-        nDVs = len(xDV[designVarName])
+        nDVs = 0
+        parameters = None
+        if self.DVGeo is not None:
+            xDV = self.DVGeo.getValues()
+            if designVarName in xDV:
+                nDVs = len(xDV[designVarName])
+                parameters = xDV[designVarName].copy(order="C")
+
+        if designVarName in self.internalDV:
+            nDVs = len(self.internalDV[designVarName]["init"])
+            parameters = self.internalDV[designVarName]["value"]
 
         nParameters = self.solver.getNRegressionParameters()
         if nDVs != nParameters:
@@ -2300,7 +2389,6 @@ def calcTotalDerivsRegPar(self, objFuncName, designVarName, dFScaling=1.0, accum
         xvArray = self.vec2Array(self.xvVec)
         wArray = self.vec2Array(self.wVec)
         seedArray = self.vec2Array(self.adjVectors[objFuncName])
-        parameters = xDV[designVarName].copy(order="C")
         totalDerivArray = np.zeros(nDVs)
         dFdRegPar = np.zeros(nDVs)
 
@@ -2371,6 +2459,9 @@ def solveAdjointUnsteady(self):
         # self.solverAD once to get the proper oldTime level for unsteady adjoint
         self.solverAD.calcPrimalResidualStatistics("calc".encode())
 
+        # call the internal design var function to update DASolver parameters
+        self.runInternalDVFunc()
+
         # calc the total number of time instances
         # we assume the adjoint is for deltaT to endTime
         # but users can also prescribed a custom time range
@@ -2590,6 +2681,9 @@ def solveAdjoint(self):
             if self.getOption("useAD")["mode"] == "reverse":
                 self.solverAD.updateOFField(self.wVec)
 
+        # call the internal design var function to update DASolver parameters
+        self.runInternalDVFunc()
+
         self.adjointFail = 0
 
         # calculate dRdWT

src/adjoint/DARegression/DARegression.C

@@ -41,6 +41,8 @@ DARegression::DARegression(
 
     regSubDict.readEntry<label>("printInputRange", printInputRange_);
 
+    regSubDict.readEntry<scalar>("defaultOutputValue", defaultOutputValue_);
+
     active_ = regSubDict.getLabel("active");
 
     if (modelType_ == "neuralNetwork")
@@ -355,12 +357,12 @@ void DARegression::checkOutput(volScalarField& outputField)
     {
         if (std::isnan(outputField[cellI]))
         {
-            outputField[cellI] = 0;
+            outputField[cellI] = defaultOutputValue_;
             isNaN = 1;
         }
         if (std::isinf(outputField[cellI]))
         {
-            outputField[cellI] = 0;
+            outputField[cellI] = defaultOutputValue_;
             isInf = 1;
         }
         if (outputField[cellI] > outputUpperBound_)
@@ -376,17 +378,17 @@ void DARegression::checkOutput(volScalarField& outputField)
     }
     if (isBounded == 1)
     {
-        Info << "************* Warning! output values are bounded. ******************" << endl;
+        Info << "************* Warning! output values are bounded between " << outputLowerBound_ << " and " << outputUpperBound_ << endl;
     }
 
     if (isNaN == 1)
     {
-        Info << "************* Warning! output values have nan and were set to zeros ******************" << endl;
+        Info << "************* Warning! output values have nan and are set to " << defaultOutputValue_ << endl;
     }
 
     if (isInf == 1)
     {
-        Info << "************* Warning! output values have inf and were set to zeros ******************" << endl;
+        Info << "************* Warning! output values have inf and are set to " << defaultOutputValue_ << endl;
     }
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/adjoint/DARegression/DARegression.H

@@ -90,6 +90,9 @@ protected:
     /// whether to print the input range info this is used to scale the input
     label printInputRange_;
 
+    /// default output values
+    scalar defaultOutputValue_;
+
 public:
     /// Constructors
     DARegression(

tests/runTests_DAPimpleFoam.py

@@ -226,7 +226,7 @@ def regModel(val, geo):
 nParameters = DASolver.solver.getNRegressionParameters()
 
 parameter0 = np.ones(nParameters) * 0.1
-DVGeo.addGlobalDV("parameter", parameter0, regModel, lower=-100.0, upper=100.0, scale=1.0)
+DASolver.addInternalDV("parameter", parameter0, regModel, lower=-100.0, upper=100.0, scale=1.0)
 
 DASolver.setDVGeo(DVGeo)
 mesh = USMesh(options=meshOptions, comm=gcomm)
@@ -257,10 +257,12 @@ def regModel(val, geo):
 else:
     DASolver.runColoring()
     xDV = DVGeo.getValues()
+    iDV = DASolver.getInternalDVDict()
+    allDV = {**xDV, **iDV}
     funcs = {}
-    funcs, fail = optFuncs.calcObjFuncValues(xDV)
+    funcs, fail = optFuncs.calcObjFuncValues(allDV)
     funcsSens = {}
-    funcsSens, fail = optFuncs.calcObjFuncSens(xDV, funcs)
+    funcsSens, fail = optFuncs.calcObjFuncSens(allDV, funcs)
 
     parameterNormU = np.linalg.norm(funcsSens["CD"]["parameter"])
     funcsSens["CD"]["parameter"] = parameterNormU

tests/runTests_DASimpleFoam.py

@@ -162,7 +162,7 @@
 nTwists = DVGeo.addRefAxis("bodyAxis", xFraction=0.25, alignIndex="k")
 
 
-def alpha(val, geo):
+def alpha(val, DASolver):
     aoa = val[0] * np.pi / 180.0
     inletU = [float(U0 * np.cos(aoa)), float(U0 * np.sin(aoa)), 0]
     DASolver.setOption("primalBC", {"U0": {"variable": "U", "patches": ["inout"], "value": inletU}})
@@ -212,7 +212,6 @@ def actuator(val, geo):
 indexList = pts[1:4, 1, 0].flatten()
 PS = geo_utils.PointSelect("list", indexList)
 DVGeo.addLocalDV("shapey", lower=-1.0, upper=1.0, axis="y", scale=1.0, pointSelect=PS)
-DVGeo.addGlobalDV("alpha", [alpha0], alpha, lower=-10.0, upper=10.0, scale=1.0)
 # actuator
 DVGeo.addGlobalDV(
     "actuator",
@@ -225,6 +224,7 @@ def actuator(val, geo):
 
 # DAFoam
 DASolver = PYDAFOAM(options=aeroOptions, comm=gcomm)
+DASolver.addInternalDV("alpha", [alpha0], alpha, lower=-10.0, upper=10.0, scale=1.0)
 DASolver.setDVGeo(DVGeo)
 mesh = USMesh(options=meshOptions, comm=gcomm)
 DASolver.printFamilyList()
@@ -253,16 +253,18 @@ def actuator(val, geo):
 else:
     DASolver.runColoring()
     xDV = DVGeo.getValues()
+    iDV = DASolver.getInternalDVDict()
+    allDV = {**xDV, **iDV}
     funcs = {}
-    funcs, fail = optFuncs.calcObjFuncValues(xDV)
+    funcs, fail = optFuncs.calcObjFuncValues(allDV)
     # test getForces
     forces = DASolver.getForces()
     fNorm = np.linalg.norm(forces.flatten())
     fNormSum = gcomm.allreduce(fNorm, op=MPI.SUM)
     funcs["forces"] = fNormSum
 
     funcsSens = {}
-    funcsSens, fail = optFuncs.calcObjFuncSens(xDV, funcs)
+    funcsSens, fail = optFuncs.calcObjFuncSens(allDV, funcs)
     if gcomm.rank == 0:
         reg_write_dict(funcs, 1e-8, 1e-10)
         reg_write_dict(funcsSens, 1e-4, 1e-6)