Added a regression model class (#547)

* Added DARegression and its derivs.

* Called correctBC for regModel.

* Fixed a missing init bug.

* Added an option to not time vol for the volSum obj.

* Changed the DAReg init.

* Fixed a bug.

* Added tests and fixed an issue in compiling.

* Fixed an issuse for Reg.

dafoam/pyDAFoam.py

@@ -508,6 +508,17 @@ def __init__(self):
         ## constrainHbyA back to the primal and adjoint solvers.
         self.useConstrainHbyA = False
 
+        ## parameters for regression models
+        self.regressionModel = {
+            "active": False,
+            "modelType": "neuralNetwork",
+            "inputNames": ["None"],
+            "outputName": "None",
+            "hiddenLayerNeurons": [0],
+            "outputShift": 0.0,
+            "outputScale": 1.0,
+        }
+
         # *********************************************************************************************
         # ************************************ Advance Options ****************************************
         # *********************************************************************************************
@@ -2264,6 +2275,44 @@ def calcTotalDerivsACT(self, objFuncName, designVarName, designVarType, dFScalin
             else:
                 self.adjTotalDeriv[objFuncName][designVarName][i] = totalDerivSeq[i]
 
+    def calcTotalDerivsRegPar(self, objFuncName, designVarName, accumulateTotal=False):
+
+        xDV = self.DVGeo.getValues()
+        nDVs = len(xDV[designVarName])
+
+        nParameters = self.solver.getNRegressionParameters()
+        if nDVs != nParameters:
+            raise Error("number of parameters not valid!")
+
+        # We assume  dFdRegPar is always zero
+        # call the total deriv
+
+        xvArray = self.vec2Array(self.xvVec)
+        wArray = self.vec2Array(self.wVec)
+        seedArray = self.vec2Array(self.adjVectors[objFuncName])
+        parameters = xDV[designVarName].copy(order="C")
+        productArray = np.zeros(nDVs)
+
+        # calculate dRdFieldT*Psi and save it to totalDeriv
+        self.solverAD.calcdRdRegParTPsiAD(xvArray, wArray, parameters, seedArray, productArray)
+        # all reduce because parameters is a global DV
+        productArray = self.comm.allreduce(productArray, op=MPI.SUM)
+
+        # totalDeriv = dFdRegPar - dRdRegParT*psi
+        productArray *= -1.0
+
+        # assign the total derivative to self.adjTotalDeriv
+        if self.adjTotalDeriv[objFuncName][designVarName] is None:
+            self.adjTotalDeriv[objFuncName][designVarName] = np.zeros(nDVs, self.dtype)
+
+        # NOTE: productArray is already in Seq
+
+        for i in range(nDVs):
+            if accumulateTotal is True:
+                self.adjTotalDeriv[objFuncName][designVarName][i] += productArray[i]
+            else:
+                self.adjTotalDeriv[objFuncName][designVarName][i] = productArray[i]
+
     def solveAdjointUnsteady(self):
         """
         Run adjoint solver to compute the total derivs for unsteady solvers
@@ -2274,6 +2323,8 @@ def solveAdjointUnsteady(self):
 
         wVec: vector that contains all the state variables
 
+        designVariable: a dictionary that contain the design variable values
+
         Output:
         -------
         self.adjTotalDeriv: the dict contains all the total derivative vectors
@@ -2462,6 +2513,8 @@ def solveAdjointUnsteady(self):
                         elif designVarDict[designVarName]["designVarType"] == "Field":
                             fieldType = designVarDict[designVarName]["fieldType"]
                             self.calcTotalDerivsField(objFuncName, designVarName, fieldType, dFScaling, True)
+                        elif designVarDict[designVarName]["designVarType"] == "RegPar":
+                            self.calcTotalDerivsRegPar(objFuncName, designVarName, True)
                         else:
                             raise Error("designVarType not valid!")
 
@@ -2751,6 +2804,15 @@ def solveAdjoint(self):
                             self.calcTotalDerivsField(objFuncName, designVarName, fieldType)
                 else:
                     raise Error("For Field design variable type, we only support useAD->mode=reverse")
+            ################### RegPar: regression model's parameters ###################
+            elif designVarDict[designVarName]["designVarType"] == "RegPar":
+                if self.getOption("useAD")["mode"] == "reverse":
+                    # loop over all objectives
+                    for objFuncName in objFuncDict:
+                        if objFuncName in self.objFuncNames4Adj:
+                            self.calcTotalDerivsRegPar(objFuncName, designVarName)
+                else:
+                    raise Error("For Field design variable type, we only support useAD->mode=reverse")
             else:
                 raise Error("designVarType %s not supported!" % designVarDict[designVarName]["designVarType"])
 
@@ -3151,7 +3213,7 @@ def _initSolver(self):
         adjMode = self.getOption("unsteadyAdjoint")["mode"]
         if adjMode == "hybrid":
             self.initTimeInstanceMats()
-        
+
         Info("Init solver done! ElapsedClockTime %f s" % self.solver.getElapsedClockTime())
         Info("Init solver done!. ElapsedCpuTime %f s" % self.solver.getElapsedCpuTime())
 
@@ -4002,6 +4064,23 @@ def _initOption(self, name, value):
                 "Expected data type is %-47s \n "
                 "Received data type is %-47s" % (name, self.defaultOptions[name][0], type(value))
             )
+
+    def setRegressionParameter(self, idx, val):
+        """
+        Update the regression parameters
+
+        Parameters
+        ----------
+        idx: int
+            the index of the parameter
+
+        val, float
+            the parameter value to set
+        """
+
+        self.solver.setRegressionParameter(idx, val)
+        if self.getOption("useAD")["mode"] in ["forward", "reverse"]:
+            self.solverAD.setRegressionParameter(idx, val)
 
     def setFieldValue4GlobalCellI(self, fieldName, val, globalCellI, compI=0):
         """

src/adjoint/DAModel/DAModel.C

@@ -292,6 +292,28 @@ void DAModel::getTurbProdTerm(scalarList& prodTerm) const
 
 }
 
+void DAModel::getTurbProdOverDestruct(scalarList& PoD) const
+{
+    /*
+    Description: 
+        Return the value of the production/destruction term from the turbulence model 
+    */
+
+#ifndef SolidDASolver
+    if (hasTurbulenceModel_)
+    {
+        DATurbulenceModel& daTurb = const_cast<DATurbulenceModel&>(
+            mesh_.thisDb().lookupObject<DATurbulenceModel>("DATurbulenceModel"));
+        daTurb.getTurbProdOverDestruct(PoD);
+    }
+
+    if (hasRadiationModel_)
+    {
+    }
+#endif
+
+}
+
 #ifdef CompressibleFlow
 const fluidThermo& DAModel::getThermo() const
 {

src/adjoint/DAModel/DAModel.H

@@ -105,9 +105,12 @@ public:
     /// update intermediate variables that are dependent on the model states
     void updateIntermediateVariables();
 
-    /// return the value of the production term from the turbulence model 
+    /// return the value of the production term from the turbulence model
     void getTurbProdTerm(scalarList& prodTerm) const;
 
+    /// return the value of the destruction term from the turbulence model
+    void getTurbProdOverDestruct(scalarList& PoD) const;
+
 #ifndef SolidDASolver
     /// get a reference to DATurbulenceModel
     const DATurbulenceModel& getDATurbulenceModel() const;
@@ -117,7 +120,6 @@ public:
     /// return thermo_ from DATurbulenceModel
     const fluidThermo& getThermo() const;
 #endif
-
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.C

@@ -518,6 +518,27 @@ void DASpalartAllmaras::getFvMatrixFields(
     upper = nuTildaEqn.upper();
     lower = nuTildaEqn.lower();
 }
+
+void DASpalartAllmaras::getTurbProdOverDestruct(scalarList& PoD) const
+{
+    /*
+    Description:
+        Return the value of the production over destruction term from the turbulence model 
+    */
+
+    const volScalarField chi(this->chi());
+    const volScalarField fv1(this->fv1(chi));
+
+    const volScalarField Stilda(this->Stilda(chi, fv1));
+
+    volScalarField P = Cb1_ * phase_ * rho_ * Stilda * nuTilda_;
+    volScalarField D = Cw1_ * phase_ * rho_ * fw(Stilda) * sqr(nuTilda_ / y_);
+
+    forAll(P, cellI)
+    {
+        PoD[cellI] = P[cellI] / D[cellI];
+    }
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.H

@@ -146,6 +146,9 @@ public:
         scalarField& diag,
         scalarField& upper,
         scalarField& lower);
+
+    /// return the value of the destruction term from the turbulence model
+    virtual void getTurbProdOverDestruct(scalarList& PoD) const;
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.C

@@ -737,13 +737,36 @@ void DASpalartAllmarasFv3::getFvMatrixFields(
             - Cb2_ / sigmaNut_ * phase_ * rho_ * magSqr(fvc::grad(nuTilda_))
         == Cb1_ * phase_ * rho_ * Stilda * nuTilda_ * betaFI_
             - fvm::Sp(Cw1_ * phase_ * rho_ * fw(Stilda) * nuTilda_ / sqr(y_), nuTilda_));
-    
+
     nuTildaEqn.relax();
 
     diag = nuTildaEqn.D();
     upper = nuTildaEqn.upper();
     lower = nuTildaEqn.lower();
 }
+
+void DASpalartAllmarasFv3::getTurbProdOverDestruct(scalarList& PoD) const
+{
+    /*
+    Description:
+        Return the value of the production over destruction term from the turbulence model 
+    */
+
+    const volScalarField chi(this->chi());
+    const volScalarField fv1(this->fv1(chi));
+
+    const volScalarField Stilda(
+        this->fv3(chi, fv1) * ::sqrt(2.0) * mag(skew(fvc::grad(U_)))
+        + this->fv2(chi, fv1) * nuTilda_ / sqr(kappa_ * y_));
+
+    volScalarField P = Cb1_ * phase_ * rho_ * Stilda * nuTilda_;
+    volScalarField D = Cw1_ * phase_ * rho_ * fw(Stilda) * sqr(nuTilda_ / y_);
+
+    forAll(P, cellI)
+    {
+        PoD[cellI] = P[cellI] / D[cellI];
+    }
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.H

@@ -168,6 +168,9 @@ public:
         scalarField& diag,
         scalarField& upper,
         scalarField& lower);
+
+    /// return the value of the destruction term from the turbulence model
+    virtual void getTurbProdOverDestruct(scalarList& PoD) const;
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.C

@@ -487,7 +487,19 @@ void DATurbulenceModel::getTurbProdTerm(scalarList& prodTerm) const
         Return the value of the production term from the turbulence model 
     */
 
-    FatalErrorIn("DATurbulenceModel::getSAProdTerm")
+    FatalErrorIn("DATurbulenceModel::getTurbProdTerm")
+        << "Child class not implemented!"
+        << abort(FatalError);
+}
+
+void DATurbulenceModel::getTurbProdOverDestruct(scalarList& PoD) const
+{
+    /*
+    Description:
+        Return the value of the production over destruction term from the turbulence model 
+    */
+
+    FatalErrorIn("DATurbulenceModel::getTurbProdOverDestruct")
         << "Child class not implemented!"
         << abort(FatalError);
 }

src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.H

@@ -211,9 +211,12 @@ public:
     /// solve the residual equations and update the state
     virtual void correct(label printToScreen) = 0;
 
-    /// return the value of the production term from the Spalart Allmaras model
+    /// return the value of the production term from the turbulence model
     virtual void getTurbProdTerm(scalarList& prodTerm) const;
 
+    /// return the value of the destruction term from the turbulence model
+    virtual void getTurbProdOverDestruct(scalarList& PoD) const;
+
     /// dev terms
     tmp<volSymmTensorField> devRhoReff() const;
 

src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C

@@ -51,6 +51,8 @@ DAObjFuncVariableVolSum::DAObjFuncVariableVolSum(
 
     objFuncDict_.readEntry<label>("divByTotalVol", divByTotalVol_);
 
+    multiplyVol_ = objFuncDict_.lookupOrDefault<label>("multiplyVol", 1);
+
     if (daIndex.adjStateNames.found(varName_))
     {
         objFuncConInfo_ = {{varName_}};
@@ -114,7 +116,11 @@ void DAObjFuncVariableVolSum::calcObjFunc(
         forAll(objFuncCellSources, idxI)
         {
             const label& cellI = objFuncCellSources[idxI];
-            scalar volume = mesh_.V()[cellI];
+            scalar volume = 1.0;
+            if (multiplyVol_)
+            {
+                volume = mesh_.V()[cellI];
+            }
             if (isSquare_)
             {
                 objFuncCellValues[idxI] = scale_ * volume * var[cellI] * var[cellI];
@@ -133,7 +139,11 @@ void DAObjFuncVariableVolSum::calcObjFunc(
         forAll(objFuncCellSources, idxI)
         {
             const label& cellI = objFuncCellSources[idxI];
-            scalar volume = mesh_.V()[cellI];
+            scalar volume = 1.0;
+            if (multiplyVol_)
+            {
+                volume = mesh_.V()[cellI];
+            }
             if (isSquare_)
             {
                 objFuncCellValues[idxI] = scale_ * volume * var[cellI][component_] * var[cellI][component_];

src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.H

@@ -44,6 +44,9 @@ protected:
     /// whether to normalize the volume sum with the total volume
     label divByTotalVol_;
 
+    /// whether to multiply the variable by the volume
+    label multiplyVol_;
+
 public:
     TypeName("variableVolSum");
     // Constructors