diff --git a/src/adjoint/DAFunction/DAFunctionVariableVolSum.C b/src/adjoint/DAFunction/DAFunctionVariableVolSum.C
index c920e375..55d1359f 100755
--- a/src/adjoint/DAFunction/DAFunctionVariableVolSum.C
+++ b/src/adjoint/DAFunction/DAFunctionVariableVolSum.C
@@ -41,6 +41,10 @@ DAFunctionVariableVolSum::DAFunctionVariableVolSum(
     multiplyVol_ = functionDict_.lookupOrDefault<label>("multiplyVol", 1);
 
     divByTotalVol_ = functionDict_.lookupOrDefault<label>("divByTotalVol", 0);
+
+    invertField_ = functionDict_.lookupOrDefault<label>("invertField", 0);
+
+    invertVal_ = functionDict_.lookupOrDefault<label>("invertVal", 1); // Binary inversion from 0 to 1, vice versa
 }
 
 /// calculate the value of objective function
@@ -57,7 +61,7 @@ scalar DAFunctionVariableVolSum::calcFunction()
 
     const objectRegistry& db = mesh_.thisDb();
 
-    scalar totalVol = 1.0;
+    scalar totalVol = 0.0;
 
     if (divByTotalVol_)
     {
@@ -67,6 +71,10 @@ scalar DAFunctionVariableVolSum::calcFunction()
         }
         reduce(totalVol, sumOp<scalar>());
     }
+    else
+    {
+        totalVol = 1.0;
+    }
 
     if (varType_ == "scalar")
     {
@@ -75,40 +83,52 @@ scalar DAFunctionVariableVolSum::calcFunction()
         forAll(cellSources_, idxI)
         {
             const label& cellI = cellSources_[idxI];
+            scalar value = var[cellI];
             scalar volume = 1.0;
             if (multiplyVol_)
             {
                 volume = mesh_.V()[cellI];
             }
+            if (invertField_)
+            {
+
+                value = invertVal_ - value;
+            }
             if (isSquare_)
             {
-                functionValue += scale_ * volume * var[cellI] * var[cellI];
+                functionValue += scale_ * volume * value * value;
             }
             else
             {
-                functionValue += scale_ * volume * var[cellI];
+                functionValue += scale_ * volume * value;
             }
         }
     }
     else if (varType_ == "vector")
     {
+       
         const volVectorField& var = db.lookupObject<volVectorField>(varName_);
         // calculate mass
         forAll(cellSources_, idxI)
         {
             const label& cellI = cellSources_[idxI];
+            scalar value = var[cellI][index_];
             scalar volume = 1.0;
             if (multiplyVol_)
             {
                 volume = mesh_.V()[cellI];
             }
+            if (invertField_)
+            {
+                value = invertVal_ - value;
+            }
             if (isSquare_)
             {
-                functionValue += scale_ * volume * var[cellI][index_] * var[cellI][index_];
+                functionValue += scale_ * volume * value * value;
             }
             else
             {
-                functionValue += scale_ * volume * var[cellI][index_];
+                functionValue += scale_ * volume * value;
             }
         }
     }
@@ -122,8 +142,14 @@ scalar DAFunctionVariableVolSum::calcFunction()
     // need to reduce the sum of force across all processors
     reduce(functionValue, sumOp<scalar>());
 
+    Info << "Raw functionValue = " << functionValue << endl;
+    Info << "Total Volume      = " << totalVol << endl;
+    Info << "Ratio             = " << functionValue/totalVol << endl;
+    Info << "volume            = " << totalVol << endl;
+
     functionValue /= totalVol;
 
+
     // check if we need to calculate refDiff.
     this->calcRefVar(functionValue);
 
diff --git a/src/adjoint/DAFunction/DAFunctionVariableVolSum.H b/src/adjoint/DAFunction/DAFunctionVariableVolSum.H
old mode 100644
new mode 100755
index 551ea2bc..02115ead
--- a/src/adjoint/DAFunction/DAFunctionVariableVolSum.H
+++ b/src/adjoint/DAFunction/DAFunctionVariableVolSum.H
@@ -46,6 +46,12 @@ protected:
     /// whether to multiply the variable by the volume
     label multiplyVol_;
 
+    /// whether to invert the binary scalar field
+    label invertField_;
+
+    /// scalar value to invert the field by
+    scalar invertVal_;
+
 public:
     TypeName("variableVolSum");
     // Constructors
diff --git a/src/adjoint/DAFunction/DAFunctionViscousDissipation.C b/src/adjoint/DAFunction/DAFunctionViscousDissipation.C
new file mode 100644
index 00000000..7eb61ac5
--- /dev/null
+++ b/src/adjoint/DAFunction/DAFunctionViscousDissipation.C
@@ -0,0 +1,114 @@
+#include "DAFunctionViscousDissipation.H"
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(DAFunctionViscousDissipation, 0);
+addToRunTimeSelectionTable(
+    DAFunction,
+    DAFunctionViscousDissipation,
+    dictionary);
+
+DAFunctionViscousDissipation::DAFunctionViscousDissipation(
+    const fvMesh& mesh,
+    const DAOption& daOption,
+    const DAModel& daModel,
+    const DAIndex& daIndex,
+    const word functionName)
+    : DAFunction(
+        mesh,
+        daOption,
+        daModel,
+        daIndex,
+        functionName),
+      daTurb_(daModel.getDATurbulenceModel())
+{
+}
+
+scalar DAFunctionViscousDissipation::calcFunction()
+{
+    /*
+    Description:
+        Calculate the total power dissipation (Reference: https://doi.org/10.1002/nme.1468)
+
+            J = \int [ 2*muEff*D:D + alphaPorosity*(U·U) ] dV
+
+        where
+
+            D = symm(grad(U))
+
+        The first term represents viscous dissipation while the second
+        represents Brinkman (Darcy) dissipation.
+
+    Output:
+        functionValue
+    */
+
+    scalar functionValue = 0.0;
+
+    const objectRegistry& db = mesh_.thisDb();
+
+    //-------------------------------------------------------------
+    // Flow variables
+    //-------------------------------------------------------------
+
+    const volVectorField& U =
+        db.lookupObject<volVectorField>("U");
+
+    const volScalarField& alpha =
+        db.lookupObject<volScalarField>("alpha");
+
+    //-------------------------------------------------------------
+    // Velocity gradient
+    //-------------------------------------------------------------
+
+    tmp<volTensorField> tGradU = fvc::grad(U);
+
+    //-------------------------------------------------------------
+    // Symmetric strain-rate tensor
+    //-------------------------------------------------------------
+
+    volSymmTensorField D(
+        symm(tGradU()));
+
+    //-------------------------------------------------------------
+    // Effective viscosity
+    //-------------------------------------------------------------
+
+    tmp<volScalarField> tNuEff = daTurb_.nuEff();
+
+    //-------------------------------------------------------------
+    // Integrate objective
+    //-------------------------------------------------------------
+
+    forAll(cellSources_, idxI)
+    {
+        const label cellI = cellSources_[idxI];
+
+        scalar viscousTerm =
+            2.0
+            * tNuEff()[cellI]
+            * (D[cellI] && D[cellI]);
+
+        scalar brinkmanTerm =
+            alpha[cellI]
+            * magSqr(U[cellI]);
+
+        functionValue +=
+            scale_
+            * (viscousTerm + brinkmanTerm)
+            * mesh_.V()[cellI];
+    }
+
+    //-------------------------------------------------------------
+    // Parallel reduction
+    //-------------------------------------------------------------
+
+    reduce(functionValue, sumOp<scalar>());
+
+    this->calcRefVar(functionValue);
+
+    return functionValue;
+}
+
+}
diff --git a/src/adjoint/DAFunction/DAFunctionViscousDissipation.H b/src/adjoint/DAFunction/DAFunctionViscousDissipation.H
new file mode 100644
index 00000000..ee8dbdb6
--- /dev/null
+++ b/src/adjoint/DAFunction/DAFunctionViscousDissipation.H
@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+
+    DAFoam : Discrete Adjoint with OpenFOAM
+
+    Description:
+        Power (viscous) dissipation objective function
+    Usage:
+        "function": {
+        "viscousDissipation": {
+            "type": "DAViscousDissipation",
+            "scale": 1.0,
+            "source": "allCells"
+        }
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DAFunctionViscousDissipation_H
+#define DAFunctionViscousDissipation_H
+
+#include "DAFunction.H"
+#include "addToRunTimeSelectionTable.H"
+
+namespace Foam
+{
+
+class DAFunctionViscousDissipation
+    : public DAFunction
+{
+
+protected:
+    /// Turbulence model
+    const DATurbulenceModel& daTurb_;
+
+public:
+    TypeName("DAViscousDissipation");
+
+    DAFunctionViscousDissipation(
+        const fvMesh& mesh,
+        const DAOption& daOption,
+        const DAModel& daModel,
+        const DAIndex& daIndex,
+        const word functionName);
+
+    virtual ~DAFunctionViscousDissipation()
+    {
+    }
+
+    virtual scalar calcFunction();
+};
+
+}
+
+#endif
diff --git a/src/adjoint/Make/files b/src/adjoint/Make/files
index 87385b23..b148264f 100755
--- a/src/adjoint/Make/files
+++ b/src/adjoint/Make/files
@@ -43,6 +43,7 @@ DAFunction/DAFunctionMeshQualityKS.C
 DAFunction/DAFunctionVonMisesStressKS.C
 DAFunction/DAFunctionResidualNorm.C
 DAFunction/DAFunctionFieldMax.C
+DAFunction/DAFunctionViscousDissipation.C
 
 DAFvSource/DAFvSource.C
 DAFvSource/DAFvSourceActuatorDisk.C