Added an option to use time-dependent ref for refDiffSquare.

Author: friedenhe

src/adjoint/DAObjFunc/DAObjFunc.C

@@ -59,7 +59,7 @@ DAObjFunc::DAObjFunc(
     calcRefDiffSquare_ = objFuncDict_.lookupOrDefault<label>("calcRefDiffSquare", 0);
     if (calcRefDiffSquare_)
     {
-        objFuncDict_.readEntry<scalar>("ref", ref_);
+        objFuncDict_.readEntry<scalarList>("ref", ref_);
     }
 }
 
@@ -318,6 +318,29 @@ scalar DAObjFunc::getObjFuncValue()
     return objFuncValue_;
 }
 
+void DAObjFunc::calcRefDiffSquare(scalar& objFuncValue)
+{
+    /*
+    Description:
+        Call the variable difference with respect to a given reference and take a square of it.
+        This can be used in FIML. This function is for calcRefDiffSquare_ == 1
+    */
+
+    if (calcRefDiffSquare_)
+    {
+        if (ref_.size() == 1)
+        {
+            objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
+        }
+        else
+        {
+            label idxI = runTime.timeIndex() - 1;
+            objFuncValue = (objFuncValue - ref_[idxI]) * (objFuncValue - ref_[idxI]);
+        }
+    }
+
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/adjoint/DAObjFunc/DAObjFunc.H

@@ -95,9 +95,8 @@ protected:
     /// whether to calculate (obj-ref)^2
     label calcRefDiffSquare_;
 
-    /// if calcRefDiffSquare_ is True, set the reference value
-    scalar ref_;
-
+    /// if calcRefDiffSquare_ is True, set the reference value list
+    scalarList ref_;
 
 public:
     /// Runtime type information
@@ -243,6 +242,9 @@ public:
 
     /// whether to compute reference coefficients for special objFunc treatment such as totalPressureRatio
     label calcRefCoeffs = 1;
+
+    /// calculate the varDiffSquare
+    void calcRefDiffSquare(scalar& objFuncValue);
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/adjoint/DAObjFunc/DAObjFuncCenterOfPressure.C

@@ -128,10 +128,8 @@ void DAObjFuncCenterOfPressure::calcObjFunc(
 
     objFuncValue = scale_ * (weightedPressure / totalPressure + (center_ & axis_));
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncForce.C

@@ -191,10 +191,8 @@ void DAObjFuncForce::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncLocation.C

@@ -219,10 +219,8 @@ void DAObjFuncLocation::calcObjFunc(
                                           << abort(FatalError);
     }
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncMass.C

@@ -96,10 +96,8 @@ void DAObjFuncMass::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncMassFlowRate.C

@@ -119,10 +119,8 @@ void DAObjFuncMassFlowRate::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncMeshQualityKS.C

@@ -219,10 +219,8 @@ void DAObjFuncMeshQualityKS::calcObjFunc(
 
     objFuncValue = log(objFuncValue) / coeffKS_;
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncMoment.C

@@ -165,10 +165,8 @@ void DAObjFuncMoment::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncPatchMean.C

@@ -132,10 +132,8 @@ void DAObjFuncPatchMean::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncPower.C

@@ -168,10 +168,8 @@ void DAObjFuncPower::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncTotalPressure.C

@@ -132,10 +132,8 @@ void DAObjFuncTotalPressure::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncTotalPressureRatio.C

@@ -217,10 +217,8 @@ void DAObjFuncTotalPressureRatio::calcObjFunc(
         }
     }
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncTotalTemperatureRatio.C

@@ -207,10 +207,8 @@ void DAObjFuncTotalTemperatureRatio::calcObjFunc(
         }
     }
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C

@@ -167,10 +167,8 @@ void DAObjFuncVariableVolSum::calcObjFunc(
 
     objFuncValue /= totalVol;
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncVonMisesStressKS.C

@@ -121,10 +121,8 @@ void DAObjFuncVonMisesStressKS::calcObjFunc(
 
     objFuncValue = log(objValTmp) / coeffKS_;
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

src/adjoint/DAObjFunc/DAObjFuncWallHeatFlux.C

@@ -246,10 +246,8 @@ void DAObjFuncWallHeatFlux::calcObjFunc(
     // need to reduce the sum of force across all processors
     reduce(objFuncValue, sumOp<scalar>());
 
-    if (calcRefDiffSquare_)
-    {
-        objFuncValue = (objFuncValue - ref_) * (objFuncValue - ref_);
-    }
+    // check if we need to calculate refDiffSquare.
+    this->calcRefDiffSquare(objFuncValue);
 
     return;
 }

tests/runTests_DASimpleFoamAD.py

@@ -131,7 +131,7 @@
                 "scale": 1.0 / (0.5 * U0 * U0 * A0 * 1.0),
                 "addToAdjoint": True,
                 "calcRefDiffSquare": True,
-                "ref": 0.1
+                "ref": [0.1]
             }
         },
     },