Added an option to compute objective function difference square. (#588)

* Updated the version.

* Added an option to consider refDiffSquare for all objs.

* Fixed the test and a potential issue in dFdAOA-analytical.

* Fixed a test script and ref.

* Updated test ref.

dafoam/pyDAFoam.py

@@ -11,7 +11,7 @@
 
 """
 
-__version__ = "3.1.0"
+__version__ = "3.1.1"
 
 import subprocess
 import os
@@ -3031,6 +3031,23 @@ def calcdFdAOAAnalytical(self, objFuncName, dFdAOA):
 
         # if it is a forceObj, use the analytical approach to calculate dFdAOA, otherwise set it to zero
         if isForceObj == 1:
+            # First, we need to check if a pair of lift and drag forces
+            # is defined. If not, return an error.
+            nParallelToFlows = 0
+            nNormalToFlows = 0
+            for objFuncName in objFuncDict:
+                for objFuncPart in objFuncDict[objFuncName]:
+                    if objFuncDict[objFuncName][objFuncPart]["type"] == "force":
+                        if objFuncDict[objFuncName][objFuncPart]["directionMode"] == "parallelToFlow":
+                            nParallelToFlows += 1
+                        if objFuncDict[objFuncName][objFuncPart]["directionMode"] == "normalToFlow":
+                            nNormalToFlows += 1
+
+            if nParallelToFlows != 1 or nNormalToFlows != 1:
+                raise Error(
+                    "calcdFdAOAAnalytical supports only one pair of force: one with normalToFlow and the other with parallelToFlow"
+                )
+
             # loop over all objectives again to find the neededMode
             # Note that if the neededMode == "parallelToFlow", we need to add a minus sign
             for objFuncNameNeeded in objFuncDict:

src/adjoint/DAObjFunc/DAObjFunc.C

@@ -55,6 +55,12 @@ DAObjFunc::DAObjFunc(
     {
         objFuncCellValues_[idxI] = 0.0;
     }
+
+    calcRefDiffSquare_ = objFuncDict_.lookupOrDefault<label>("calcRefDiffSquare", 0);
+    if (calcRefDiffSquare_)
+    {
+        objFuncDict_.readEntry<scalar>("ref", ref_);
+    }
 }
 
 // * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //

src/adjoint/DAObjFunc/DAObjFunc.H

@@ -92,6 +92,13 @@ protected:
     /// the connectivity information for the objective function
     List<List<word>> objFuncConInfo_;
 
+    /// whether to calculate (obj-ref)^2
+    label calcRefDiffSquare_;
+
+    /// if calcRefDiffSquare_ is True, set the reference value
+    scalar ref_;
+
+
 public:
     /// Runtime type information
     TypeName("DAObjFunc");

src/adjoint/DAObjFunc/DAObjFuncCenterOfPressure.C

@@ -128,6 +128,11 @@ void DAObjFuncCenterOfPressure::calcObjFunc(
 
     objFuncValue = scale_ * (weightedPressure / totalPressure + (center_ & axis_));
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/adjoint/DAObjFunc/DAObjFuncForce.C

@@ -191,6 +191,11 @@ void DAObjFuncForce::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 

src/adjoint/DAObjFunc/DAObjFuncLocation.C

@@ -219,6 +219,11 @@ void DAObjFuncLocation::calcObjFunc(
                                           << abort(FatalError);
     }
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 

src/adjoint/DAObjFunc/DAObjFuncMass.C

@@ -96,6 +96,11 @@ void DAObjFuncMass::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 

src/adjoint/DAObjFunc/DAObjFuncMassFlowRate.C

@@ -119,6 +119,11 @@ void DAObjFuncMassFlowRate::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 

src/adjoint/DAObjFunc/DAObjFuncMeshQualityKS.C

@@ -219,6 +219,11 @@ void DAObjFuncMeshQualityKS::calcObjFunc(
 
     objFuncValue = log(objFuncValue) / coeffKS_;
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 

src/adjoint/DAObjFunc/DAObjFuncMoment.C

@@ -165,6 +165,11 @@ void DAObjFuncMoment::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/adjoint/DAObjFunc/DAObjFuncPatchMean.C

@@ -132,6 +132,11 @@ void DAObjFuncPatchMean::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 

src/adjoint/DAObjFunc/DAObjFuncPower.C

@@ -168,6 +168,11 @@ void DAObjFuncPower::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/adjoint/DAObjFunc/DAObjFuncTotalPressure.C

@@ -132,6 +132,11 @@ void DAObjFuncTotalPressure::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 

src/adjoint/DAObjFunc/DAObjFuncTotalPressureRatio.C

@@ -217,6 +217,11 @@ void DAObjFuncTotalPressureRatio::calcObjFunc(
         }
     }
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 

src/adjoint/DAObjFunc/DAObjFuncTotalTemperatureRatio.C

@@ -207,6 +207,11 @@ void DAObjFuncTotalTemperatureRatio::calcObjFunc(
         }
     }
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 

src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C

@@ -167,6 +167,11 @@ void DAObjFuncVariableVolSum::calcObjFunc(
 
     objFuncValue /= totalVol;
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 

src/adjoint/DAObjFunc/DAObjFuncVonMisesStressKS.C

@@ -121,6 +121,11 @@ void DAObjFuncVonMisesStressKS::calcObjFunc(
 
     objFuncValue = log(objValTmp) / coeffKS_;
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 

src/adjoint/DAObjFunc/DAObjFuncWallHeatFlux.C

@@ -246,6 +246,11 @@ void DAObjFuncWallHeatFlux::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
+    if (calcRefDiffSquare_)
+    {
+        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+    }
+
     return;
 }
 

tests/refs/DAFoam_Test_DASimpleFoamADRef.txt

@@ -2,6 +2,8 @@ Dictionary Key: CD
 @value  0.02043001230171681 1e-08 1e-10
 Dictionary Key: CL
 @value   0.4874216250833033 1e-08 1e-10
+Dictionary Key: CMZDIFF
+@value   0.0115816476130927 1e-08 1e-10
 Dictionary Key: COP
 @value   0.2298218575123644 1e-08 1e-10
 Dictionary Key: R1
@@ -30,6 +32,15 @@ Dictionary Key: shapey
 @value    1.418731716621286 1e-05 1e-07
 Dictionary Key: twist
 @value  0.09008667655578639 1e-05 1e-07
+Dictionary Key: CMZDIFF
+Dictionary Key: alpha
+@value 0.0005090113344677443 1e-05 1e-07
+Dictionary Key: shapey
+@value -0.001151982015258751 1e-05 1e-07
+@value -0.03924921702616808 1e-05 1e-07
+@value  -0.1132153003029267 1e-05 1e-07
+Dictionary Key: twist
+@value 0.0004899695826971005 1e-05 1e-07
 Dictionary Key: COP
 Dictionary Key: alpha
 @value -0.002070288134429892 1e-05 1e-07

tests/refs/DAFoam_Test_MphysAeroFieldRef.txt

@@ -4,6 +4,6 @@ Dictionary Key: CL
 Dictionary Key: alpha
 @value    0.435012393483863 0.0001 1e-06
 Dictionary Key: aoa
-@value  0.09767145516644697 0.0001 1e-06
+@value  0.0973750318897765 0.0001 1e-06
 Dictionary Key: twist
 @value  0.09790920083964022 0.0001 1e-06

tests/runTests_DASimpleFoamAD.py

@@ -110,6 +110,30 @@
                 "addToAdjoint": True,
             }
         },
+        "CMZ": {
+            "part1": {
+                "type": "moment",
+                "source": "patchToFace",
+                "patches": ["wing"],
+                "axis": [0.0, 0.0, 1.0],
+                "center": [0.25, 0.0, 0.05],
+                "scale": 1.0 / (0.5 * U0 * U0 * A0 * 1.0),
+                "addToAdjoint": False
+            }
+        },
+        "CMZDIFF": {
+            "part1": {
+                "type": "moment",
+                "source": "patchToFace",
+                "patches": ["wing"],
+                "axis": [0.0, 0.0, 1.0],
+                "center": [0.25, 0.0, 0.05],
+                "scale": 1.0 / (0.5 * U0 * U0 * A0 * 1.0),
+                "addToAdjoint": True,
+                "calcRefDiffSquare": True,
+                "ref": 0.1
+            }
+        },
     },
     "normalizeStates": {"U": U0, "p": U0 * U0 / 2.0, "k": k0, "omega": omega0, "phi": 1.0},
     "adjPartDerivFDStep": {"State": 1e-6, "FFD": 1e-3, "ACTD": 1.0e-3},

tests/runTests_MphysAeroField.py

@@ -43,6 +43,17 @@
         "transport:nu": 1.5e-5,
     },
     "objFunc": {
+        "CD": {
+            "part1": {
+                "type": "force",
+                "source": "patchToFace",
+                "patches": ["wing"],
+                "directionMode": "parallelToFlow",
+                "alphaName": "aoa",
+                "scale": 1.0 / (0.5 * U0 * U0 * A0),
+                "addToAdjoint": True,
+            }
+        },
         "CL": {
             "part1": {
                 "type": "force",