From a941eaf42d632fba4084f550f94a10a8ce99d4ab Mon Sep 17 00:00:00 2001
From: Ping He <phe@iastate.edu>
Date: Thu, 27 Jun 2024 18:02:36 -0500
Subject: [PATCH] Allowed heat source to move with derivs (#656)

* Added heat source parameters as DVs and location snapCenter option.

* Fixed some bugs and added the derivs for the heatSource pars.

* Fixed a bug.
---
 dafoam/mphys/mphys_dafoam.py                  |  50 +++
 dafoam/pyDAFoam.py                            |   2 +-
 src/adjoint/DAFvSource/DAFvSource.C           |  96 +++---
 src/adjoint/DAFvSource/DAFvSource.H           |  22 +-
 .../DAFvSource/DAFvSourceActuatorDisk.C       |  38 +--
 src/adjoint/DAFvSource/DAFvSourceHeatSource.C | 147 +++++++--
 src/adjoint/DAFvSource/DAFvSourceHeatSource.H |   5 +-
 src/adjoint/DAObjFunc/DAObjFuncLocation.C     |  36 +-
 src/adjoint/DAObjFunc/DAObjFuncLocation.H     |   6 +
 .../DAResidual/DAResidualHeatTransferFoam.C   |  13 +
 .../DAResidual/DAResidualHeatTransferFoam.H   |   2 +
 .../DAHeatTransferFoam/DAHeatTransferFoam.C   |   4 -
 src/adjoint/DASolver/DASolver.C               | 310 ++++++++++++++++++
 src/adjoint/DASolver/DASolver.H               |  15 +
 src/pyDASolvers/DASolvers.H                   |  22 ++
 src/pyDASolvers/pyDASolvers.pyx               |   8 +
 .../DAFoam_Test_DAHeatTransferFoamRef.txt     |  12 +-
 .../refs/DAFoam_Test_MphysAeroThermalRef.txt  |  12 +
 tests/runTests_DAHeatTransferFoam.py          |  49 ++-
 tests/runTests_MphysAeroThermal.py            |  37 ++-
 20 files changed, 774 insertions(+), 112 deletions(-)

diff --git a/dafoam/mphys/mphys_dafoam.py b/dafoam/mphys/mphys_dafoam.py
index 35c44978..15522ed2 100644
--- a/dafoam/mphys/mphys_dafoam.py
+++ b/dafoam/mphys/mphys_dafoam.py
@@ -657,6 +657,11 @@ def setup(self):
                 if "comps" in list(designVariables[dvName].keys()):
                     nACTDVars = len(designVariables[dvName]["comps"])
                 self.add_input(dvName, distributed=False, shape=nACTDVars, tags=["mphys_coupling"])
+            elif dvType == "HSC":  # add heat source parameter variables
+                nHSCVars = 9
+                if "comps" in list(designVariables[dvName].keys()):
+                    nHSCVars = len(designVariables[dvName]["comps"])
+                self.add_input(dvName, distributed=False, shape=nHSCVars, tags=["mphys_coupling"])
             elif dvType == "Field":  # add field variables
                 # user can prescribe whether the field var is distributed. Default is True
                 distributed = True
@@ -888,6 +893,26 @@ def apply_linear(self, inputs, outputs, d_inputs, d_outputs, d_residuals, mode):
                             d_inputs[inputName] += ACTDBarSub
                         else:
                             d_inputs[inputName] += ACTDBar
+                    
+                    # compute [dRdHSC]^T*Psi using reverse mode AD
+                    elif self.dvType[inputName] == "HSC":
+                        prodVec = PETSc.Vec().create(self.comm)
+                        prodVec.setSizes((PETSc.DECIDE, 9), bsize=1)
+                        prodVec.setFromOptions()
+                        DASolver.solverAD.calcdRdHSCTPsiAD(
+                            DASolver.xvVec, DASolver.wVec, resBarVec, inputName.encode(), prodVec
+                        )
+                        # we will convert the MPI prodVec to seq array for all procs
+                        HSCBar = DASolver.convertMPIVec2SeqArray(prodVec)
+                        if "comps" in list(designVariables[inputName].keys()):
+                            nHSCVars = len(designVariables[inputName]["comps"])
+                            HSCBarSub = np.zeros(nHSCVars, "d")
+                            for i in range(nHSCVars):
+                                comp = designVariables[inputName]["comps"][i]
+                                HSCBarSub[i] = HSCBar[comp]
+                            d_inputs[inputName] += HSCBarSub
+                        else:
+                            d_inputs[inputName] += HSCBar
 
                     # compute dRdFieldT*Psi using reverse mode AD
                     elif self.dvType[inputName] == "Field":
@@ -1251,6 +1276,11 @@ def setup(self):
                 if "comps" in list(designVariables[dvName].keys()):
                     nACTDVars = len(designVariables[dvName]["comps"])
                 self.add_input(dvName, distributed=False, shape=nACTDVars, tags=["mphys_coupling"])
+            elif dvType == "HSC":  # add heat source parameter variables
+                nHSCVars = 9
+                if "comps" in list(designVariables[dvName].keys()):
+                    nHSCVars = len(designVariables[dvName]["comps"])
+                self.add_input(dvName, distributed=False, shape=nHSCVars, tags=["mphys_coupling"])
             elif dvType == "Field":  # add field variables
                 # user can prescribe whether the field var is distributed. Default is True
                 distributed = True
@@ -1447,6 +1477,26 @@ def compute_jacvec_product(self, inputs, d_inputs, d_outputs, mode):
                             d_inputs[inputName] += ACTDBarSub * fBar
                         else:
                             d_inputs[inputName] += ACTDBar * fBar
+                    
+                    # compute dFdHSC
+                    elif self.dvType[inputName] == "HSC":
+                        dFdHSC = PETSc.Vec().create(self.comm)
+                        dFdHSC.setSizes((PETSc.DECIDE, 9), bsize=1)
+                        dFdHSC.setFromOptions()
+                        DASolver.solverAD.calcdFdHSCAD(
+                            DASolver.xvVec, DASolver.wVec, objFuncName.encode(), inputName.encode(), dFdHSC
+                        )
+                        # we will convert the MPI dFdHSC to seq array for all procs
+                        HSCBar = DASolver.convertMPIVec2SeqArray(dFdHSC)
+                        if "comps" in list(designVariables[inputName].keys()):
+                            nHSCVars = len(designVariables[inputName]["comps"])
+                            HSCBarSub = np.zeros(nHSCVars, "d")
+                            for i in range(nHSCVars):
+                                comp = designVariables[inputName]["comps"][i]
+                                HSCBarSub[i] = HSCBar[comp]
+                            d_inputs[inputName] += HSCBarSub * fBar
+                        else:
+                            d_inputs[inputName] += HSCBar * fBar
 
                     # compute dFdField
                     elif self.dvType[inputName] == "Field":
diff --git a/dafoam/pyDAFoam.py b/dafoam/pyDAFoam.py
index 0d66923f..84b337f0 100755
--- a/dafoam/pyDAFoam.py
+++ b/dafoam/pyDAFoam.py
@@ -11,7 +11,7 @@
 
 """
 
-__version__ = "3.1.2"
+__version__ = "3.1.3"
 
 import subprocess
 import os
diff --git a/src/adjoint/DAFvSource/DAFvSource.C b/src/adjoint/DAFvSource/DAFvSource.C
index a436ebe0..01ea675f 100755
--- a/src/adjoint/DAFvSource/DAFvSource.C
+++ b/src/adjoint/DAFvSource/DAFvSource.C
@@ -135,49 +135,67 @@ void DAFvSource::syncDAOptionToActuatorDVs()
 
     // now we need to initialize actuatorDiskDVs_
     dictionary fvSourceSubDict = daOption_.getAllOptions().subDict("fvSource");
-    word diskName0 = fvSourceSubDict.toc()[0];
 
-    word type0 = fvSourceSubDict.subDict(diskName0).getWord("type");
-
-    if (type0 == "actuatorDisk")
+    forAll(fvSourceSubDict.toc(), idxI)
     {
-        word source0 = fvSourceSubDict.subDict(diskName0).getWord("source");
+        word diskName = fvSourceSubDict.toc()[idxI];
+        // sub dictionary with all parameters for this disk
+        dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
+        word type = diskSubDict.getWord("type");
+        word source = diskSubDict.getWord("source");
+
+        if (type == "actuatorDisk" && source == "cylinderAnnulusSmooth")
+        {
 
-        if (source0 == "cylinderAnnulusSmooth")
+            // now read in all parameters for this actuator disk
+            scalarList centerList;
+            diskSubDict.readEntry<scalarList>("center", centerList);
+
+            scalarList dirList;
+            diskSubDict.readEntry<scalarList>("direction", dirList);
+
+            // we have 13 design variables for each disk
+            scalarList dvList(13);
+            dvList[0] = centerList[0];
+            dvList[1] = centerList[1];
+            dvList[2] = centerList[2];
+            dvList[3] = dirList[0];
+            dvList[4] = dirList[1];
+            dvList[5] = dirList[2];
+            dvList[6] = diskSubDict.getScalar("innerRadius");
+            dvList[7] = diskSubDict.getScalar("outerRadius");
+            dvList[8] = diskSubDict.getScalar("scale");
+            dvList[9] = diskSubDict.getScalar("POD");
+            dvList[10] = diskSubDict.getScalar("expM");
+            dvList[11] = diskSubDict.getScalar("expN");
+            dvList[12] = diskSubDict.getScalar("targetThrust");
+
+            // set actuatorDiskDVs_
+            actuatorDiskDVs_.set(diskName, dvList);
+        }
+        else if (type == "heatSource" && source == "cylinderSmooth")
         {
-            forAll(fvSourceSubDict.toc(), idxI)
-            {
-                word diskName = fvSourceSubDict.toc()[idxI];
-
-                // sub dictionary with all parameters for this disk
-                dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
-
-                // now read in all parameters for this actuator disk
-                scalarList centerList;
-                diskSubDict.readEntry<scalarList>("center", centerList);
-
-                scalarList dirList;
-                diskSubDict.readEntry<scalarList>("direction", dirList);
-
-                // we have 13 design variables for each disk
-                scalarList dvList(13);
-                dvList[0] = centerList[0];
-                dvList[1] = centerList[1];
-                dvList[2] = centerList[2];
-                dvList[3] = dirList[0];
-                dvList[4] = dirList[1];
-                dvList[5] = dirList[2];
-                dvList[6] = diskSubDict.getScalar("innerRadius");
-                dvList[7] = diskSubDict.getScalar("outerRadius");
-                dvList[8] = diskSubDict.getScalar("scale");
-                dvList[9] = diskSubDict.getScalar("POD");
-                dvList[10] = diskSubDict.getScalar("expM");
-                dvList[11] = diskSubDict.getScalar("expN");
-                dvList[12] = diskSubDict.getScalar("targetThrust");
-
-                // set actuatorDiskDVs_
-                actuatorDiskDVs_.set(diskName, dvList);
-            }
+            // now read in all parameters for this actuator disk
+            scalarList centerList;
+            diskSubDict.readEntry<scalarList>("center", centerList);
+
+            scalarList axisList;
+            diskSubDict.readEntry<scalarList>("axis", axisList);
+
+            // we have 13 design variables for each disk
+            scalarList dvList(9);
+            dvList[0] = centerList[0];
+            dvList[1] = centerList[1];
+            dvList[2] = centerList[2];
+            dvList[3] = axisList[0];
+            dvList[4] = axisList[1];
+            dvList[5] = axisList[2];
+            dvList[6] = diskSubDict.getScalar("radius");
+            dvList[7] = diskSubDict.getScalar("length");
+            dvList[8] = diskSubDict.getScalar("power");
+
+            // set actuatorDiskDVs_
+            actuatorDiskDVs_.set(diskName, dvList);
         }
     }
 }
diff --git a/src/adjoint/DAFvSource/DAFvSource.H b/src/adjoint/DAFvSource/DAFvSource.H
index e693f125..296577c0 100755
--- a/src/adjoint/DAFvSource/DAFvSource.H
+++ b/src/adjoint/DAFvSource/DAFvSource.H
@@ -128,9 +128,25 @@ public:
     /// calculate fvSource based on the latest actuatorDVs
     void updateFvSource()
     {
-        volVectorField& fvSource = const_cast<volVectorField&>(
-            mesh_.thisDb().lookupObject<volVectorField>("fvSource"));
-        this->calcFvSource(fvSource);
+        if (mesh_.thisDb().foundObject<volVectorField>("fvSource"))
+        {
+            volVectorField& fvSource = const_cast<volVectorField&>(
+                mesh_.thisDb().lookupObject<volVectorField>("fvSource"));
+
+            this->calcFvSource(fvSource);
+        }
+        else if (mesh_.thisDb().foundObject<volScalarField>("fvSource"))
+        {
+            volScalarField& fvSource = const_cast<volScalarField&>(
+                mesh_.thisDb().lookupObject<volScalarField>("fvSource"));
+
+            this->calcFvSource(fvSource);
+        }
+        else
+        {
+            FatalErrorIn("") << "fvSource not found! "
+                             << abort(FatalError);
+        }
     }
 
     /// synchronize the values in DAOption and actuatorDiskDVs_
diff --git a/src/adjoint/DAFvSource/DAFvSourceActuatorDisk.C b/src/adjoint/DAFvSource/DAFvSourceActuatorDisk.C
index 4b9f43a2..b950667f 100755
--- a/src/adjoint/DAFvSource/DAFvSourceActuatorDisk.C
+++ b/src/adjoint/DAFvSource/DAFvSourceActuatorDisk.C
@@ -103,21 +103,15 @@ void DAFvSourceActuatorDisk::calcFvSource(volVectorField& fvSource)
 
     dictionary fvSourceSubDict = allOptions.subDict("fvSource");
 
-    word diskName0 = fvSourceSubDict.toc()[0];
-    word source0 = fvSourceSubDict.subDict(diskName0).getWord("source");
-
-    if (source0 == "cylinderAnnulusToCell")
+    forAll(fvSourceSubDict.toc(), idxI)
     {
+        word diskName = fvSourceSubDict.toc()[idxI];
+        dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
+        word source = diskSubDict.getWord("source");
 
-        // loop over all the cell indices for all actuator disks
-        forAll(fvSourceCellIndices_.toc(), idxI)
+        if (source == "cylinderAnnulusToCell")
         {
-
-            // name of this disk
-            word diskName = fvSourceCellIndices_.toc()[idxI];
-
-            // sub dictionary with all parameters for this disk
-            dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
+            // loop over all the cell indices for all actuator disks
 
             // now read in all parameters for this actuator disk
             scalarList point1;
@@ -213,14 +207,8 @@ void DAFvSourceActuatorDisk::calcFvSource(volVectorField& fvSource)
                 }
             }
         }
-    }
-    else if (source0 == "cylinderAnnulusSmooth")
-    {
-
-        forAll(fvSourceSubDict.toc(), idxI)
+        else if (source == "cylinderAnnulusSmooth")
         {
-            word diskName = fvSourceSubDict.toc()[idxI];
-            dictionary diskSubDict = fvSourceSubDict.subDict(diskName);
 
             vector center = {actuatorDiskDVs_[diskName][0], actuatorDiskDVs_[diskName][1], actuatorDiskDVs_[diskName][2]};
             vector dirNorm = {actuatorDiskDVs_[diskName][3], actuatorDiskDVs_[diskName][4], actuatorDiskDVs_[diskName][5]};
@@ -432,12 +420,12 @@ void DAFvSourceActuatorDisk::calcFvSource(volVectorField& fvSource)
                 }
             }
         }
-    }
-    else
-    {
-        FatalErrorIn("calcFvSourceCells") << "source: " << source0 << " not supported!"
-                                          << "Options are: cylinderAnnulusToCell and cylinderAnnulusSmooth!"
-                                          << abort(FatalError);
+        else
+        {
+            FatalErrorIn("calcFvSourceCells") << "source: " << source << " not supported!"
+                                              << "Options are: cylinderAnnulusToCell and cylinderAnnulusSmooth!"
+                                              << abort(FatalError);
+        }
     }
 
     fvSource.correctBoundaryConditions();
diff --git a/src/adjoint/DAFvSource/DAFvSourceHeatSource.C b/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
index 3fff3688..69ba2abf 100755
--- a/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
+++ b/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
@@ -24,6 +24,8 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
     const DAIndex& daIndex)
     : DAFvSource(modelType, mesh, daOption, daModel, daIndex)
 {
+    printInterval_ = daOption.getOption<label>("printInterval");
+
     const dictionary& allOptions = daOption_.getAllOptions();
 
     dictionary fvSourceSubDict = allOptions.subDict("fvSource");
@@ -34,17 +36,18 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
         dictionary sourceSubDict = fvSourceSubDict.subDict(sourceName);
         word sourceType = sourceSubDict.getWord("source");
 
-        scalarList p1List;
-        sourceSubDict.readEntry<scalarList>("p1", p1List);
-        cylinderP1_.set(sourceName, p1List);
-        scalarList p2List;
-        sourceSubDict.readEntry<scalarList>("p2", p2List);
-        cylinderP2_.set(sourceName, p2List);
-        cylinderRadius_.set(sourceName, sourceSubDict.getScalar("radius"));
-        power_.set(sourceName, sourceSubDict.getScalar("power"));
-
         if (sourceType == "cylinderToCell")
         {
+            scalarList p1List;
+            sourceSubDict.readEntry<scalarList>("p1", p1List);
+            cylinderP1_.set(sourceName, p1List);
+            scalarList p2List;
+            sourceSubDict.readEntry<scalarList>("p2", p2List);
+            cylinderP2_.set(sourceName, p2List);
+
+            cylinderRadius_.set(sourceName, sourceSubDict.getScalar("radius"));
+            power_.set(sourceName, sourceSubDict.getScalar("power"));
+
             fvSourceCellIndices_.set(sourceName, {});
             // all available source type are in src/meshTools/sets/cellSources
             // Example of IO parameters os in applications/utilities/mesh/manipulation/topoSet
@@ -92,17 +95,26 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
                 Info << "fvSourceCellIndices " << fvSourceCellIndices_ << endl;
             }
         }
-        else if (sourceType == "cylinder")
+        else if (sourceType == "cylinderSmooth")
         {
-            // no need to do special things
+
+            // eps is a smoothing parameter, it should be the local mesh cell size in meters
+            // near the cylinder region
+            cylinderEps_.set(sourceName, sourceSubDict.getScalar("eps"));
         }
         else
         {
             FatalErrorIn("DAFvSourceHeatSource") << "source: " << sourceType << " not supported!"
-                                                 << "Options are: cylinderAnnulusToCell and cylinder!"
+                                                 << "Options are: cylinderAnnulusToCell and cylinderSmooth!"
                                                  << abort(FatalError);
         }
     }
+
+    // now we need to initialize actuatorDiskDVs_ by synchronizing the values
+    // defined in fvSource from DAOption to actuatorDiskDVs_
+    // NOTE: we need to call this function whenever we change the actuator
+    // design variables during optimization
+    this->syncDAOptionToActuatorDVs();
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -156,10 +168,11 @@ void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)
         word sourceType = sourceSubDict.getWord("source");
         // NOTE: here power should be in W. We will evenly divide the power by the
         // total volume of the source
-        scalar power = power_[sourceName];
         if (sourceType == "cylinderToCell")
         {
 
+            scalar power = power_[sourceName];
+
             // loop over all cell indices for this source and assign the source term
             // ----- first loop, calculate the total volume
             scalar totalV = 0.0;
@@ -184,20 +197,74 @@ void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)
 
             reduce(sourceTotal, sumOp<scalar>());
 
-            Info << "Total volume for " << sourceName << ": " << totalV << " m^3" << endl;
-            Info << "Total heat source for " << sourceName << ": " << sourceTotal << " W." << endl;
+            if (daOption_.getOption<word>("runStatus") == "solvePrimal")
+            {
+                if (mesh_.time().timeIndex() % printInterval_ == 0 || mesh_.time().timeIndex() == 1)
+                {
+                    Info << "Total volume for " << sourceName << ": " << totalV << " m^3" << endl;
+                    Info << "Total heat source for " << sourceName << ": " << sourceTotal << " W." << endl;
+                }
+            }
         }
-        else if (sourceType == "cylinder")
+        else if (sourceType == "cylinderSmooth")
         {
-            vector p1 = {cylinderP1_[sourceName][0], cylinderP1_[sourceName][1], cylinderP1_[sourceName][2]};
-            vector p2 = {cylinderP2_[sourceName][0], cylinderP2_[sourceName][1], cylinderP2_[sourceName][2]};
-            vector cylinderCenter = (p1 + p2) / 2.0;
-            vector cylinderDir = p2 - p1;
-            scalar cylinderLen = mag(cylinderDir);
-            vector cylinderDirNorm = cylinderDir / cylinderLen;
-            scalar radius = cylinderRadius_[sourceName];
-            scalar totalVolume = constant::mathematical::pi * radius * radius * cylinderLen;
+            vector cylinderCenter =
+                {actuatorDiskDVs_[sourceName][0], actuatorDiskDVs_[sourceName][1], actuatorDiskDVs_[sourceName][2]};
+            vector cylinderDir =
+                {actuatorDiskDVs_[sourceName][3], actuatorDiskDVs_[sourceName][4], actuatorDiskDVs_[sourceName][5]};
+            scalar radius = actuatorDiskDVs_[sourceName][6];
+            scalar cylinderLen = actuatorDiskDVs_[sourceName][7];
+            scalar power = actuatorDiskDVs_[sourceName][8];
+            vector cylinderDirNorm = cylinderDir / mag(cylinderDir);
+            scalar eps = cylinderEps_[sourceName];
+
+            // first loop, calculate the total volume
+            scalar volumeTotal = 0.0;
+            forAll(mesh_.cells(), cellI)
+            {
+                // the cell center coordinates of this cellI
+                vector cellC = mesh_.C()[cellI];
+                // cell center to disk center vector
+                vector cellC2AVec = cellC - cylinderCenter;
+                // tmp tensor for calculating the axial/radial components of cellC2AVec
+                tensor cellC2AVecE(tensor::zero);
+                cellC2AVecE.xx() = cellC2AVec.x();
+                cellC2AVecE.yy() = cellC2AVec.y();
+                cellC2AVecE.zz() = cellC2AVec.z();
+
+                // now we need to decompose cellC2AVec into axial and radial components
+                // the axial component of cellC2AVec vector
+                vector cellC2AVecA = cellC2AVecE & cylinderDirNorm;
+                // the radial component of cellC2AVec vector
+                vector cellC2AVecR = cellC2AVec - cellC2AVecA;
 
+                // the magnitude of radial component of cellC2AVecR
+                scalar cellC2AVecRLen = mag(cellC2AVecR);
+                // the magnitude of axial component of cellC2AVecR
+                scalar cellC2AVecALen = mag(cellC2AVecA);
+
+                scalar axialSmoothC = 1.0;
+                scalar radialSmoothC = 1.0;
+
+                scalar halfL = cylinderLen / 2.0;
+
+                if (cellC2AVecALen > halfL)
+                {
+                    scalar d2 = (cellC2AVecALen - halfL) * (cellC2AVecALen - halfL);
+                    axialSmoothC = exp(-d2 / eps / eps);
+                }
+
+                if (cellC2AVecRLen > radius)
+                {
+                    scalar d2 = (cellC2AVecRLen - radius) * (cellC2AVecRLen - radius);
+                    radialSmoothC = exp(-d2 / eps / eps);
+                }
+                volumeTotal += axialSmoothC * radialSmoothC * mesh_.V()[cellI];
+            }
+            reduce(volumeTotal, sumOp<scalar>());
+
+            // second loop, calculate fvSource
+            scalar sourceTotal = 0.0;
             forAll(mesh_.cells(), cellI)
             {
                 // the cell center coordinates of this cellI
@@ -221,16 +288,42 @@ void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)
                 // the magnitude of axial component of cellC2AVecR
                 scalar cellC2AVecALen = mag(cellC2AVecA);
 
-                if (cellC2AVecRLen <= radius && cellC2AVecALen <= cylinderLen / 2.0)
+                scalar axialSmoothC = 1.0;
+                scalar radialSmoothC = 1.0;
+
+                scalar halfL = cylinderLen / 2.0;
+
+                if (cellC2AVecALen > halfL)
+                {
+                    scalar d2 = (cellC2AVecALen - halfL) * (cellC2AVecALen - halfL);
+                    axialSmoothC = exp(-d2 / eps / eps);
+                }
+
+                if (cellC2AVecRLen > radius)
+                {
+                    scalar d2 = (cellC2AVecRLen - radius) * (cellC2AVecRLen - radius);
+                    radialSmoothC = exp(-d2 / eps / eps);
+                }
+
+                fvSource[cellI] += power / volumeTotal * axialSmoothC * radialSmoothC;
+                sourceTotal += power / volumeTotal * axialSmoothC * radialSmoothC * mesh_.V()[cellI];
+            }
+
+            reduce(sourceTotal, sumOp<scalar>());
+
+            if (daOption_.getOption<word>("runStatus") == "solvePrimal")
+            {
+                if (mesh_.time().timeIndex() % printInterval_ == 0 || mesh_.time().timeIndex() == 1)
                 {
-                    fvSource[cellI] += power / totalVolume;
+                    Info << "Total volume for " << sourceName << ": " << volumeTotal << " m^3" << endl;
+                    Info << "Total heat source for " << sourceName << ": " << sourceTotal << " W." << endl;
                 }
             }
         }
         else
         {
             FatalErrorIn("DAFvSourceHeatSource") << "source: " << sourceType << " not supported!"
-                                                 << "Options are: cylinderAnnulusToCell and cylinder!"
+                                                 << "Options are: cylinderAnnulusToCell and cylinderSmooth!"
                                                  << abort(FatalError);
         }
     }
diff --git a/src/adjoint/DAFvSource/DAFvSourceHeatSource.H b/src/adjoint/DAFvSource/DAFvSourceHeatSource.H
index 236106cd..4531c3c7 100755
--- a/src/adjoint/DAFvSource/DAFvSourceHeatSource.H
+++ b/src/adjoint/DAFvSource/DAFvSourceHeatSource.H
@@ -35,9 +35,10 @@ protected:
     HashTable<scalarList> cylinderP2_;
     HashTable<scalar> cylinderRadius_;
     HashTable<scalar> power_;
+    HashTable<scalar> cylinderEps_;
 
-    /// calculate DAFvSourceHeatSource::fvSourceCellIndices_
-    void calcFvSourceCellIndices(const word sourceName, HashTable<labelList>& fvSourceCellIndices);
+    /// print interval for primal and adjoint
+    label printInterval_;
 
 public:
     TypeName("heatSource");
diff --git a/src/adjoint/DAObjFunc/DAObjFuncLocation.C b/src/adjoint/DAObjFunc/DAObjFuncLocation.C
index 1b1e4bcd..8f2245b5 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncLocation.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncLocation.C
@@ -65,6 +65,17 @@ DAObjFuncLocation::DAObjFuncLocation(
         center_ = {centerRead[0], centerRead[1], centerRead[2]};
     }
 
+    snapCenter2Cell_ = objFuncDict_.lookupOrDefault<label>("snapCenter2Cell", 0);
+    if (snapCenter2Cell_)
+    {
+        point centerPoint = {center_[0], center_[1], center_[2]};
+        snappedCenterCellI_ = mesh_.findCell(centerPoint);
+        if (snappedCenterCellI_ < 0)
+        {
+            FatalErrorIn(" ") << "can not find cell for center " << centerPoint << abort(FatalError);
+        }
+    }
+
     if (mode_ == "maxRadius")
     {
         // we need to identify the patchI and faceI that has maxR
@@ -158,7 +169,13 @@ void DAObjFuncLocation::calcObjFunc(
             const label patchI = daIndex_.bFacePatchI[bFaceI];
             const label faceI = daIndex_.bFaceFaceI[bFaceI];
 
-            vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center_;
+            vector center = center_;
+            if (snapCenter2Cell_)
+            {
+                center = mesh_.C()[snappedCenterCellI_];
+            }
+
+            vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center;
 
             tensor faceCTensor(tensor::zero);
             faceCTensor.xx() = faceC.x();
@@ -200,7 +217,13 @@ void DAObjFuncLocation::calcObjFunc(
             const label patchI = daIndex_.bFacePatchI[bFaceI];
             const label faceI = daIndex_.bFaceFaceI[bFaceI];
 
-            vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center_;
+            vector center = center_;
+            if (snapCenter2Cell_)
+            {
+                center = mesh_.C()[snappedCenterCellI_];
+            }
+
+            vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center;
 
             tensor faceCTensor(tensor::zero);
             faceCTensor.xx() = faceC.x();
@@ -234,7 +257,14 @@ void DAObjFuncLocation::calcObjFunc(
         scalar radius = 0.0;
         if (maxRPatchI_ >= 0 && maxRFaceI_ >= 0)
         {
-            vector faceC = mesh_.Cf().boundaryField()[maxRPatchI_][maxRFaceI_] - center_;
+
+            vector center = center_;
+            if (snapCenter2Cell_)
+            {
+                center = mesh_.C()[snappedCenterCellI_];
+            }
+
+            vector faceC = mesh_.Cf().boundaryField()[maxRPatchI_][maxRFaceI_] - center;
 
             tensor faceCTensor(tensor::zero);
             faceCTensor.xx() = faceC.x();
diff --git a/src/adjoint/DAObjFunc/DAObjFuncLocation.H b/src/adjoint/DAObjFunc/DAObjFuncLocation.H
index 23026e7b..7dd1a140 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncLocation.H
+++ b/src/adjoint/DAObjFunc/DAObjFuncLocation.H
@@ -46,6 +46,12 @@ protected:
     /// maxRadius patchI
     label maxRFaceI_ = -1;
 
+    /// whether to snap the center to a cell in the mesh if yes the center will move with the mesh
+    label snapCenter2Cell_ = 0;
+
+    /// the cell index for the center if snapCenter2Cell_ = 1
+    label snappedCenterCellI_ = -1;
+
 public:
     TypeName("location");
     // Constructors
diff --git a/src/adjoint/DAResidual/DAResidualHeatTransferFoam.C b/src/adjoint/DAResidual/DAResidualHeatTransferFoam.C
index becf4e2a..0a628abc 100755
--- a/src/adjoint/DAResidual/DAResidualHeatTransferFoam.C
+++ b/src/adjoint/DAResidual/DAResidualHeatTransferFoam.C
@@ -42,6 +42,12 @@ DAResidualHeatTransferFoam::DAResidualHeatTransferFoam(
             "k",
             dimPower / dimLength / dimTemperature,
             transportProperties));
+    
+    const dictionary& allOptions = daOption.getAllOptions();
+    if (allOptions.subDict("fvSource").toc().size() != 0)
+    {
+        hasFvSource_ = 1;
+    }
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -72,6 +78,13 @@ void DAResidualHeatTransferFoam::calcResiduals(const dictionary& options)
         URes_, pRes_, TRes_, phiRes_: residual field variables
     */
 
+    if (hasFvSource_)
+    {
+        DAFvSource& daFvSource(const_cast<DAFvSource&>(
+            mesh_.thisDb().lookupObject<DAFvSource>("DAFvSource")));
+        daFvSource.calcFvSource(fvSource_);
+    }
+
     fvScalarMatrix TEqn(
         fvm::laplacian(kPtr_(), T_)
         + fvSource_);
diff --git a/src/adjoint/DAResidual/DAResidualHeatTransferFoam.H b/src/adjoint/DAResidual/DAResidualHeatTransferFoam.H
index 8af823af..4cea3f32 100755
--- a/src/adjoint/DAResidual/DAResidualHeatTransferFoam.H
+++ b/src/adjoint/DAResidual/DAResidualHeatTransferFoam.H
@@ -37,6 +37,8 @@ protected:
 
     autoPtr<dimensionedScalar> kPtr_;
 
+    label hasFvSource_;
+
 public:
     TypeName("DAHeatTransferFoam");
     // Constructors
diff --git a/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.C b/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.C
index 285f0eaa..bc02aaaf 100755
--- a/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.C
+++ b/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.C
@@ -57,10 +57,6 @@ void DAHeatTransferFoam::initSolver()
         word fvSourceType = allOptions.subDict("fvSource").subDict(sourceName).getWord("type");
         daFvSourcePtr_.reset(DAFvSource::New(
             fvSourceType, mesh, daOptionPtr_(), daModelPtr_(), daIndexPtr_()));
-        
-        // NOTE: we calculate the fvSource only once
-        volScalarField& fvSource = fvSourcePtr_();
-        daFvSourcePtr_->calcFvSource(fvSource);
     }
 }
 
diff --git a/src/adjoint/DASolver/DASolver.C b/src/adjoint/DASolver/DASolver.C
index ea92d912..6626f01c 100644
--- a/src/adjoint/DASolver/DASolver.C
+++ b/src/adjoint/DASolver/DASolver.C
@@ -6870,6 +6870,185 @@ void DASolver::calcdFdACTAD(
 #endif
 }
 
+void DASolver::calcdFdHSCAD(
+    const Vec xvVec,
+    const Vec wVec,
+    const word objFuncName,
+    const word designVarName,
+    Vec dFdHSC)
+{
+#ifdef CODI_AD_REVERSE
+    /*
+    Description:
+        Compute dFdHSC using reverse-mode AD
+        HSC: heat source with cylinders
+    
+    Input:
+
+        xvVec: the volume mesh coordinate vector
+
+        wVec: the state variable vector
+
+        objFuncName: the name of the objective function
+
+        designVarName: name of the design variable
+    
+    Output:
+        dFdHSC: dF/dHSC
+    */
+
+    Info << "Calculating dFdHSC using reverse-mode AD" << endl;
+
+    VecZeroEntries(dFdHSC);
+
+    label nDVs = 9;
+
+    // first check if the input is valid
+    dictionary dvSubDict = daOptionPtr_->getAllOptions().subDict("designVar").subDict(designVarName);
+    word designVarType = dvSubDict.getWord("designVarType");
+    if (designVarType == "HSC")
+    {
+        DAFvSource& fvSource = const_cast<DAFvSource&>(
+            meshPtr_->thisDb().lookupObject<DAFvSource>("DAFvSource"));
+
+        word heatSourceName = dvSubDict.getWord("heatSourceName");
+        dictionary fvSourceSubDict = daOptionPtr_->getAllOptions().subDict("fvSource");
+        word source = fvSourceSubDict.subDict(heatSourceName).getWord("source");
+        if (source == "cylinderSmooth")
+        {
+            this->updateOFField(wVec);
+            this->updateOFMesh(xvVec);
+
+            // get the subDict for this objective function
+            dictionary objFuncSubDict =
+                daOptionPtr_->getAllOptions().subDict("objFunc").subDict(objFuncName);
+
+            // loop over all parts for this objFuncName
+            forAll(objFuncSubDict.toc(), idxJ)
+            {
+                // get the subDict for this part
+                word objFuncPart = objFuncSubDict.toc()[idxJ];
+
+                // get objFunc from daObjFuncPtrList_
+                label objIndx = this->getObjFuncListIndex(objFuncName, objFuncPart);
+                DAObjFunc& daObjFunc = daObjFuncPtrList_[objIndx];
+
+                // get the design variable vals
+                scalarList actDVList(nDVs);
+                for (label i = 0; i < nDVs; i++)
+                {
+                    actDVList[i] = fvSource.getActuatorDVs(heatSourceName, i);
+                }
+
+                // reset tape
+                this->globalADTape_.reset();
+                // activate tape, start recording
+                this->globalADTape_.setActive();
+                // register  the input
+                for (label i = 0; i < nDVs; i++)
+                {
+                    this->globalADTape_.registerInput(actDVList[i]);
+                }
+                // set dv values to fvSource obj for all procs
+                for (label i = 0; i < nDVs; i++)
+                {
+                    fvSource.setActuatorDVs(heatSourceName, i, actDVList[i]);
+                }
+                // the actuatorDVs are updated, now we need to recompute fvSource
+                // this is not needed for the residual partials because fvSource
+                // will be automatically calculated in the UEqn, but for the
+                // obj partials, we need to manually recompute fvSource
+                fvSource.updateFvSource();
+
+                // update all intermediate variables and boundary conditions
+                this->updateStateBoundaryConditions();
+                // compute the objective function
+                scalar fRef = daObjFunc.getObjFuncValue();
+                // register f as the output
+                this->globalADTape_.registerOutput(fRef);
+                // stop recording
+                this->globalADTape_.setPassive();
+
+                // Note: since we used reduced objFunc, we only need to
+                // assign the seed for master proc
+                if (Pstream::master())
+                {
+                    fRef.setGradient(1.0);
+                }
+                // evaluate tape to compute derivative
+                this->globalADTape_.evaluate();
+
+                // assign the computed derivatives from the OpenFOAM variable to dFd*Part
+                Vec dFdHSCPart;
+                VecDuplicate(dFdHSC, &dFdHSCPart);
+                VecZeroEntries(dFdHSCPart);
+
+                for (label i = 0; i < nDVs; i++)
+                {
+                    PetscScalar valIn = actDVList[i].getGradient();
+                    // we need to do ADD_VALUES to get contribution from all procs
+                    VecSetValue(dFdHSCPart, i, valIn, ADD_VALUES);
+                }
+
+                VecAssemblyBegin(dFdHSCPart);
+                VecAssemblyEnd(dFdHSCPart);
+
+                // we need to add dFd*Part to dFd* because we want to sum
+                // all dFd*Part for all parts of this objFuncName.
+                VecAXPY(dFdHSC, 1.0, dFdHSCPart);
+
+                VecDestroy(&dFdHSCPart);
+
+                // need to clear adjoint and tape after the computation is done!
+                this->globalADTape_.clearAdjoints();
+                this->globalADTape_.reset();
+
+                // **********************************************************************************************
+                // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+                // **********************************************************************************************
+
+                for (label i = 0; i < nDVs; i++)
+                {
+                    this->globalADTape_.deactivateValue(actDVList[i]);
+                }
+                for (label i = 0; i < nDVs; i++)
+                {
+                    fvSource.setActuatorDVs(heatSourceName, i, actDVList[i]);
+                }
+                fvSource.updateFvSource();
+                this->updateStateBoundaryConditions();
+                fRef = daObjFunc.getObjFuncValue();
+
+                if (daOptionPtr_->getOption<label>("debug"))
+                {
+                    this->calcPrimalResidualStatistics("print");
+                    Info << objFuncName << ": " << fRef << endl;
+                }
+            }
+        }
+        else
+        {
+            FatalErrorIn("") << "source not supported. Options: cylinderAnnulusSmooth"
+                             << abort(FatalError);
+        }
+    }
+    else
+    {
+        FatalErrorIn("") << "designVarType not supported. Options: HSC"
+                         << abort(FatalError);
+    }
+
+    wordList writeJacobians;
+    daOptionPtr_->getAllOptions().readEntry<wordList>("writeJacobians", writeJacobians);
+    if (writeJacobians.found("dFdHSC") || writeJacobians.found("all"))
+    {
+        word outputName = "dFdHSC_" + objFuncName + "_" + designVarName;
+        DAUtility::writeVectorBinary(dFdHSC, outputName);
+        DAUtility::writeVectorASCII(dFdHSC, outputName);
+    }
+#endif
+}
+
 void DASolver::calcdRdActTPsiAD(
     const Vec xvVec,
     const Vec wVec,
@@ -7000,6 +7179,137 @@ void DASolver::calcdRdActTPsiAD(
 #endif
 }
 
+void DASolver::calcdRdHSCTPsiAD(
+    const Vec xvVec,
+    const Vec wVec,
+    const Vec psi,
+    const word designVarName,
+    Vec dRdHSCTPsi)
+{
+#ifdef CODI_AD_REVERSE
+    /*
+    Description:
+        Compute the matrix-vector products dRdHST^T*Psi using reverse-mode AD
+        HSC = heat source with cylinders
+    
+    Input:
+
+        xvVec: the volume mesh coordinate vector
+
+        wVec: the state variable vector
+
+        psi: the vector to multiply dRdAct
+
+        designVarName: name of the design variable
+    
+    Output:
+        dRdHSCTPsi: the matrix-vector products dRdHSC^T * Psi
+    */
+
+    Info << "Calculating [dRdHSC]^T * Psi using reverse-mode AD" << endl;
+
+    VecZeroEntries(dRdHSCTPsi);
+
+    label nDVs = 9;
+
+    dictionary dvSubDict = daOptionPtr_->getAllOptions().subDict("designVar").subDict(designVarName);
+    word designVarType = dvSubDict.getWord("designVarType");
+    if (designVarType == "HSC")
+    {
+
+        DAFvSource& fvSource = const_cast<DAFvSource&>(
+            meshPtr_->thisDb().lookupObject<DAFvSource>("DAFvSource"));
+
+        word heatSourceName = dvSubDict.getWord("heatSourceName");
+
+        dictionary fvSourceSubDict = daOptionPtr_->getAllOptions().subDict("fvSource");
+        word source = fvSourceSubDict.subDict(heatSourceName).getWord("source");
+        if (source == "cylinderSmooth")
+        {
+
+            this->updateOFField(wVec);
+            this->updateOFMesh(xvVec);
+
+            scalarList actDVList(nDVs);
+            for (label i = 0; i < nDVs; i++)
+            {
+                actDVList[i] = fvSource.getActuatorDVs(heatSourceName, i);
+            }
+
+            this->globalADTape_.reset();
+            this->globalADTape_.setActive();
+
+            for (label i = 0; i < nDVs; i++)
+            {
+                this->globalADTape_.registerInput(actDVList[i]);
+            }
+
+            // set dv values to fvSource obj for all procs
+            for (label i = 0; i < nDVs; i++)
+            {
+                fvSource.setActuatorDVs(heatSourceName, i, actDVList[i]);
+            }
+
+            // compute residuals
+            this->updateStateBoundaryConditions();
+            this->calcResiduals();
+
+            this->registerResidualOutput4AD();
+            this->globalADTape_.setPassive();
+
+            this->assignVec2ResidualGradient(psi);
+            this->globalADTape_.evaluate();
+
+            for (label i = 0; i < nDVs; i++)
+            {
+                PetscScalar valIn = actDVList[i].getGradient();
+                // we need to do ADD_VALUES to get contribution from all procs
+                VecSetValue(dRdHSCTPsi, i, valIn, ADD_VALUES);
+            }
+
+            VecAssemblyBegin(dRdHSCTPsi);
+            VecAssemblyEnd(dRdHSCTPsi);
+
+            this->globalADTape_.clearAdjoints();
+            this->globalADTape_.reset();
+
+            // **********************************************************************************************
+            // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time
+            // **********************************************************************************************
+
+            for (label i = 0; i < nDVs; i++)
+            {
+                this->globalADTape_.deactivateValue(actDVList[i]);
+            }
+            for (label i = 0; i < nDVs; i++)
+            {
+                fvSource.setActuatorDVs(heatSourceName, i, actDVList[i]);
+            }
+            this->updateStateBoundaryConditions();
+            this->calcResiduals();
+        }
+        else
+        {
+            FatalErrorIn("") << "source not supported. Options: cylinderSmooth"
+                             << abort(FatalError);
+        }
+    }
+    else
+    {
+        FatalErrorIn("") << "designVarType not supported. Options: ACTD"
+                         << abort(FatalError);
+    }
+    wordList writeJacobians;
+    daOptionPtr_->getAllOptions().readEntry<wordList>("writeJacobians", writeJacobians);
+    if (writeJacobians.found("dRdHSCTPsi") || writeJacobians.found("all"))
+    {
+        word outputName = "dRdHSCTPsi_" + designVarName;
+        DAUtility::writeVectorBinary(dRdHSCTPsi, outputName);
+        DAUtility::writeVectorASCII(dRdHSCTPsi, outputName);
+    }
+#endif
+}
+
 void DASolver::calcdForcedWAD(
     const Vec xvVec,
     const Vec wVec,
diff --git a/src/adjoint/DASolver/DASolver.H b/src/adjoint/DASolver/DASolver.H
index 864f3aea..9ac668c7 100644
--- a/src/adjoint/DASolver/DASolver.H
+++ b/src/adjoint/DASolver/DASolver.H
@@ -581,6 +581,13 @@ public:
         const Vec psi,
         const word designVarName,
         Vec dRdActTPsi);
+    
+    void calcdRdHSCTPsiAD(
+        const Vec xvVec,
+        const Vec wVec,
+        const Vec psi,
+        const word designVarName,
+        Vec dRdHSCTPsi);
 
     void calcdForcedWAD(
         const Vec xvVec,
@@ -603,6 +610,14 @@ public:
         const word objFuncName,
         const word designVarName,
         Vec dFdACT);
+    
+    /// compute dFdHSC AD
+    void calcdFdHSCAD(
+        const Vec xvVec,
+        const Vec wVec,
+        const word objFuncName,
+        const word designVarName,
+        Vec dFdHSC);
 
     /// compute dFdRegPar AD
     void calcdFdRegParAD(
diff --git a/src/pyDASolvers/DASolvers.H b/src/pyDASolvers/DASolvers.H
index 4cb4d254..0471eb52 100755
--- a/src/pyDASolvers/DASolvers.H
+++ b/src/pyDASolvers/DASolvers.H
@@ -152,6 +152,17 @@ public:
         DASolverPtr_->calcdRdActTPsiAD(xvVec, wVec, psi, designVarName, dRdActTPsi);
     }
 
+    /// compute dRdHSC^T*Psi
+    void calcdRdHSCTPsiAD(
+        const Vec xvVec,
+        const Vec wVec,
+        const Vec psi,
+        const word designVarName,
+        Vec dRdHSCTPsi)
+    {
+        DASolverPtr_->calcdRdHSCTPsiAD(xvVec, wVec, psi, designVarName, dRdHSCTPsi);
+    }
+
     /// compute dForcedW
     void calcdForcedWAD(
         const Vec xvVec,
@@ -348,6 +359,17 @@ public:
         DASolverPtr_->calcdFdACTAD(xvVec, wVec, objFuncName, designVarName, dFdACT);
     }
 
+    /// compute dFdHSC
+    void calcdFdHSCAD(
+        const Vec xvVec,
+        const Vec wVec,
+        const word objFuncName,
+        const word designVarName,
+        Vec dFdHSC)
+    {
+        DASolverPtr_->calcdFdHSCAD(xvVec, wVec, objFuncName, designVarName, dFdHSC);
+    }
+
     /// compute dFdACT
     void calcdFdACT(
         const Vec xvVec,
diff --git a/src/pyDASolvers/pyDASolvers.pyx b/src/pyDASolvers/pyDASolvers.pyx
index 7eea85df..5efc94dc 100755
--- a/src/pyDASolvers/pyDASolvers.pyx
+++ b/src/pyDASolvers/pyDASolvers.pyx
@@ -67,10 +67,12 @@ cdef extern from "DASolvers.H" namespace "Foam":
         void calcdForcedXvAD(PetscVec, PetscVec, PetscVec, PetscVec)
         void calcdAcousticsdXvAD(PetscVec, PetscVec, PetscVec, PetscVec, char*, char*)
         void calcdRdActTPsiAD(PetscVec, PetscVec, PetscVec, char*, PetscVec)
+        void calcdRdHSCTPsiAD(PetscVec, PetscVec, PetscVec, char*, PetscVec)
         void calcdForcedWAD(PetscVec, PetscVec, PetscVec, PetscVec)
         void calcdAcousticsdWAD(PetscVec, PetscVec, PetscVec, PetscVec, char*, char*)
         void calcdFdACT(PetscVec, PetscVec, char *, char*, char*, PetscVec)
         void calcdFdACTAD(PetscVec, PetscVec, char *, char*, PetscVec)
+        void calcdFdHSCAD(PetscVec, PetscVec, char *, char*, PetscVec)
         void calcdRdAOATPsiAD(PetscVec, PetscVec, PetscVec, char*, PetscVec)
         void calcdRdBCTPsiAD(PetscVec, PetscVec, PetscVec, char*, PetscVec)
         void calcdFdFFD(PetscVec, PetscVec, char *, char *, PetscVec)
@@ -263,6 +265,9 @@ cdef class pyDASolvers:
 
     def calcdRdActTPsiAD(self, Vec xvVec, Vec wVec, Vec psi, designVarName, Vec dRdActTPsi):
         self._thisptr.calcdRdActTPsiAD(xvVec.vec, wVec.vec, psi.vec, designVarName, dRdActTPsi.vec)
+    
+    def calcdRdHSCTPsiAD(self, Vec xvVec, Vec wVec, Vec psi, designVarName, Vec dRdHSCTPsi):
+        self._thisptr.calcdRdHSCTPsiAD(xvVec.vec, wVec.vec, psi.vec, designVarName, dRdHSCTPsi.vec)
 
     def calcdForcedWAD(self, Vec xvVec, Vec wVec, Vec fBarVec, Vec dForcedW):
         self._thisptr.calcdForcedWAD(xvVec.vec, wVec.vec, fBarVec.vec, dForcedW.vec)
@@ -273,6 +278,9 @@ cdef class pyDASolvers:
     def calcdFdACTAD(self, Vec xvVec, Vec wVec, objFuncName, designVarName, Vec dFdACT):
         self._thisptr.calcdFdACTAD(xvVec.vec, wVec.vec, objFuncName, designVarName, dFdACT.vec)
     
+    def calcdFdHSCAD(self, Vec xvVec, Vec wVec, objFuncName, designVarName, Vec dFdHSC):
+        self._thisptr.calcdFdHSCAD(xvVec.vec, wVec.vec, objFuncName, designVarName, dFdHSC.vec)
+    
     def calcdFdACT(self, Vec xvVec, Vec wVec, objFuncName, designVarName, designVarType, Vec dFdACT):
         self._thisptr.calcdFdACT(xvVec.vec, wVec.vec, objFuncName, designVarName, designVarType, dFdACT.vec)
     
diff --git a/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt b/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt
index e96e9016..248039f3 100755
--- a/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt
+++ b/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt
@@ -1,8 +1,14 @@
 Dictionary Key: HF_INNER
-@value    181261.4809205142 1e-08 1e-10
+@value    179999.9999999994 1e-08 1e-10
 Dictionary Key: HF_OUTER
-@value   -218738.5190794858 1e-08 1e-10
+@value   -220000.0000000005 1e-08 1e-10
+Dictionary Key: IR1KS
+@value    19.89626068314585 1e-08 1e-10
+Dictionary Key: R1
+@value    0.994617581067802 1e-08 1e-10
+Dictionary Key: R1KS
+@value    6.267981796909868 1e-08 1e-10
 Dictionary Key: TNormSum
-@value    6400011.350827799 1e-08 1e-10
+@value     6400011.35124065 1e-08 1e-10
 Dictionary Key: fail
 @value                    0 1e-08 1e-10
\ No newline at end of file
diff --git a/tests/refs/DAFoam_Test_MphysAeroThermalRef.txt b/tests/refs/DAFoam_Test_MphysAeroThermalRef.txt
index e3659018..deb2f5ce 100755
--- a/tests/refs/DAFoam_Test_MphysAeroThermalRef.txt
+++ b/tests/refs/DAFoam_Test_MphysAeroThermalRef.txt
@@ -1,4 +1,8 @@
 Dictionary Key: HFH
+Dictionary Key: power
+@value                   -0 0.0001 1e-06
+@value                   -0 0.0001 1e-06
+@value -0.003645560624152854 0.0001 1e-06
 Dictionary Key: shape
 @value   -7.000156220333093 0.0001 1e-06
 @value   -7.001072836471439 0.0001 1e-06
@@ -13,6 +17,10 @@ Dictionary Key: shape
 @value     6.49083054426032 0.0001 1e-06
 @value    6.489772648699112 0.0001 1e-06
 Dictionary Key: HF_INNER
+Dictionary Key: power
+@value                   -0 0.0001 1e-06
+@value                   -0 0.0001 1e-06
+@value   -10.00000947843464 0.0001 1e-06
 Dictionary Key: shape
 @value   -4.202340280841435 0.0001 1e-06
 @value   -4.217222098170075 0.0001 1e-06
@@ -27,6 +35,10 @@ Dictionary Key: shape
 @value    3.880687326277144 0.0001 1e-06
 @value    3.882122590063406 0.0001 1e-06
 Dictionary Key: PL
+Dictionary Key: power
+@value                   -0 0.0001 1e-06
+@value                   -0 0.0001 1e-06
+@value 2.829496338255899e-06 0.0001 1e-06
 Dictionary Key: shape
 @value     5.92445981653273 0.0001 1e-06
 @value    5.924265114818993 0.0001 1e-06
diff --git a/tests/runTests_DAHeatTransferFoam.py b/tests/runTests_DAHeatTransferFoam.py
index bc3fb2b8..3bd9727e 100755
--- a/tests/runTests_DAHeatTransferFoam.py
+++ b/tests/runTests_DAHeatTransferFoam.py
@@ -39,6 +39,45 @@
                 "addToAdjoint": True,
             }
         },
+        "R1": {
+            "part1": {
+                "type": "location",
+                "source": "patchToFace",
+                "patches": ["channel_inner"],
+                "mode": "maxRadius",
+                "axis": [0.0, 1.0, 0.0],
+                "center": [0.5, 0.05, -0.05],
+                "coeffKS": 1.0,
+                "scale": 1.0,
+                "addToAdjoint": True,
+            }
+        },
+        "R1KS": {
+            "part1": {
+                "type": "location",
+                "source": "patchToFace",
+                "patches": ["channel_inner"],
+                "mode": "maxRadiusKS",
+                "axis": [0.0, 1.0, 0.0],
+                "center": [0.5, 0.05, -0.05],
+                "coeffKS": 1.0,
+                "scale": 1.0,
+                "addToAdjoint": True,
+            }
+        },
+        "IR1KS": {
+            "part1": {
+                "type": "location",
+                "source": "patchToFace",
+                "patches": ["channel_inner"],
+                "mode": "maxInverseRadiusKS",
+                "axis": [0.0, 1.0, 0.0],
+                "center": [0.5, 0.05, -0.05],
+                "coeffKS": 1.0,
+                "scale": 1.0,
+                "addToAdjoint": True,
+            }
+        },
     },
     "debug": False,
     "primalMinResTol": 1e-12,
@@ -61,11 +100,13 @@
         },
         "source2": {
             "type": "heatSource",
-            "source": "cylinder",
-            "p1": [0.0, 0.055, 0.025],
-            "p2": [0.4, 0.055, 0.025],
+            "source": "cylinderSmooth",
+            "center": [0.2, 0.055, 0.025],
+            "axis": [1.0, 0.0, 0.0],
+            "length": 0.4,
             "radius": 0.005,
             "power": 1000.0,
+            "eps": 0.001,
         },
     },
 }
@@ -73,7 +114,7 @@
 
 DASolver()
 funcs = {}
-evalFuncs = ["HF_INNER", "HF_OUTER"]
+evalFuncs = ["HF_INNER", "HF_OUTER", "R1", "R1KS", "IR1KS"]
 DASolver.evalFunctions(funcs, evalFuncs)
 
 # test getThermal and setThermal
diff --git a/tests/runTests_MphysAeroThermal.py b/tests/runTests_MphysAeroThermal.py
index ae2248f1..48f58790 100755
--- a/tests/runTests_MphysAeroThermal.py
+++ b/tests/runTests_MphysAeroThermal.py
@@ -20,7 +20,7 @@
 
 os.chdir("./reg_test_files-main/ChannelConjugateHeat")
 if gcomm.rank == 0:
-    os.system("rm -rf 0 processor*")
+    os.system("rm -rf */processor*")
 
 # aero setup
 U0 = 10.0
@@ -178,6 +178,23 @@
     },
     "designVar": {
         "shape": {"designVarType": "FFD"},
+        "power": {
+            "designVarType": "HSC",
+            "heatSourceName": "source1",
+            "comps": [6, 7, 8],
+        },
+    },
+    "fvSource": {
+        "source1": {
+            "type": "heatSource",
+            "source": "cylinderSmooth",
+            "center": [0.2, 0.055, 0.025],
+            "axis": [1.0, 0.0, 0.0],
+            "radius": 0.005,
+            "length": 0.4,
+            "power": 10.0,
+            "eps": 0.001,
+        },
     },
 }
 
@@ -190,6 +207,15 @@
 }
 
 
+def power(val, DASolver):
+    radius = float(val[0])
+    length = float(val[1])
+    power = float(val[2])
+    # only change design variables
+    DASolver.setOption("fvSource", {"source1": {"power": power, "radius": radius, "length": length}})
+    DASolver.updateDAOption()
+
+
 class Top(Multipoint):
     def setup(self):
 
@@ -261,15 +287,21 @@ def configure(self):
         nShapes = self.geometry_aero.nom_addLocalDV(dvName="shape", axis="y", pointSelect=PS)
         nShapes = self.geometry_thermal.nom_addLocalDV(dvName="shape", axis="y", pointSelect=PS)
 
+        self.scenario.coupling.thermal.solver.add_dv_func("power", power)
+        self.scenario.thermal_post.add_dv_func("power", power)
+
         # add the design variables to the dvs component's output
         self.dvs.add_output("shape", val=np.array([0] * nShapes))
+        self.dvs.add_output("power", val=np.array([0.005, 0.4, 0.0]))
         # manually connect the dvs output to the geometry and cruise
         # sa and sst cases share the same shape
         self.connect("shape", "geometry_aero.shape")
         self.connect("shape", "geometry_thermal.shape")
+        self.connect("power", "scenario.power")
 
         # define the design variables to the top level
         self.add_design_var("shape", lower=-1, upper=1, scaler=1.0)
+        self.add_design_var("power", lower=-100.0, upper=100.0, scaler=1.0)
 
         # add objective and constraints to the top level
         self.add_objective("scenario.aero_post.HFH", scaler=1.0)
@@ -292,8 +324,11 @@ def configure(self):
     derivDict = {}
     derivDict["PL"] = {}
     derivDict["PL"]["shape"] = totals[("scenario.aero_post.functionals.PL", "dvs.shape")][0]
+    derivDict["PL"]["power"] = totals[("scenario.aero_post.functionals.PL", "dvs.power")][0]
     derivDict["HFH"] = {}
     derivDict["HFH"]["shape"] = totals[("scenario.aero_post.functionals.HFH", "dvs.shape")][0]
+    derivDict["HFH"]["power"] = totals[("scenario.aero_post.functionals.HFH", "dvs.power")][0]
     derivDict["HF_INNER"] = {}
     derivDict["HF_INNER"]["shape"] = totals[("scenario.thermal_post.functionals.HF_INNER", "dvs.shape")][0]
+    derivDict["HF_INNER"]["power"] = totals[("scenario.thermal_post.functionals.HF_INNER", "dvs.power")][0]
     reg_write_dict(derivDict, 1e-4, 1e-6)