diff --git a/src/adjoint/DAFvSource/DAFvSourceHeatSource.C b/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
index 13ec0fa8..3fff3688 100755
--- a/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
+++ b/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
@@ -24,7 +24,85 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
     const DAIndex& daIndex)
     : DAFvSource(modelType, mesh, daOption, daModel, daIndex)
 {
-    this->calcFvSourceCellIndices(fvSourceCellIndices_);
+    const dictionary& allOptions = daOption_.getAllOptions();
+
+    dictionary fvSourceSubDict = allOptions.subDict("fvSource");
+
+    forAll(fvSourceSubDict.toc(), idxI)
+    {
+        word sourceName = fvSourceSubDict.toc()[idxI];
+        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")
+        {
+            fvSourceCellIndices_.set(sourceName, {});
+            // all available source type are in src/meshTools/sets/cellSources
+            // Example of IO parameters os in applications/utilities/mesh/manipulation/topoSet
+
+            // create a topoSet
+            autoPtr<topoSet> currentSet(
+                topoSet::New(
+                    "cellSet",
+                    mesh_,
+                    "set0",
+                    IOobject::NO_READ));
+            // we need to change the min and max because they need to
+            // be of type point; however, we can't parse point type
+            // in pyDict, we need to change them here.
+
+            point p1;
+            point p2;
+            p1[0] = cylinderP1_[sourceName][0];
+            p1[1] = cylinderP1_[sourceName][1];
+            p1[2] = cylinderP1_[sourceName][2];
+            p2[0] = cylinderP2_[sourceName][0];
+            p2[1] = cylinderP2_[sourceName][1];
+            p2[2] = cylinderP2_[sourceName][2];
+
+            dictionary tmpDict;
+            tmpDict.set("p1", p1);
+            tmpDict.set("p2", p2);
+            tmpDict.set("radius", cylinderRadius_[sourceName]);
+
+            // create the source
+            autoPtr<topoSetSource> sourceSet(
+                topoSetSource::New("cylinderToCell", mesh_, tmpDict));
+
+            // add the sourceSet to topoSet
+            sourceSet().applyToSet(topoSetSource::NEW, currentSet());
+            // get the face index from currentSet, we need to use
+            // this special for loop
+            for (const label i : currentSet())
+            {
+                fvSourceCellIndices_[sourceName].append(i);
+            }
+
+            if (daOption_.getOption<label>("debug"))
+            {
+                Info << "fvSourceCellIndices " << fvSourceCellIndices_ << endl;
+            }
+        }
+        else if (sourceType == "cylinder")
+        {
+            // no need to do special things
+        }
+        else
+        {
+            FatalErrorIn("DAFvSourceHeatSource") << "source: " << sourceType << " not supported!"
+                                                 << "Options are: cylinderAnnulusToCell and cylinder!"
+                                                 << abort(FatalError);
+        }
+    }
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -43,21 +121,19 @@ void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)
                 "source1"
                 {
                     "type": "heatSource",
-                    "source": "cylinderAnnulusToCell",
+                    "source": "cylinderToCell",
                     "p1": [0.5, 0.3, 0.1], # p1 and p2 define the axis and width
                     "p2": [0.5, 0.3, 0.5], # p2-p1 should be the axis of the cylinder
-                    "innerRadius": 0.1,
-                    "outerRadius": 0.8,
+                    "radius": 0.8,
                     "power": 101.0,  # here we should prescribe the power in W
                 },
                 "source2"
                 {
                     "type": "heatSource",
-                    "source": "cylinderAnnulusToCell",
+                    "source": "cylinder",
                     "p1": [0.5, 0.3, 0.1], # p1 and p2 define the axis and width
                     "p2": [0.5, 0.3, 0.5], # p2-p1 should be the axis of the cylinder
-                    "innerRadius": 0.0,
-                    "outerRadius": 1.5,
+                    "radius": 1.5,
                     "power": 10.5, # here we should prescribe the power in W
                 }
             }
@@ -73,25 +149,17 @@ void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)
 
     dictionary fvSourceSubDict = allOptions.subDict("fvSource");
 
-    word sourceName0 = fvSourceSubDict.toc()[0];
-    word source0 = fvSourceSubDict.subDict(sourceName0).getWord("source");
-
-    if (source0 == "cylinderAnnulusToCell")
+    forAll(fvSourceSubDict.toc(), idxI)
     {
-
-        // loop over all the cell indices for all heat sources
-        forAll(fvSourceCellIndices_.toc(), idxI)
+        word sourceName = fvSourceSubDict.toc()[idxI];
+        dictionary sourceSubDict = fvSourceSubDict.subDict(sourceName);
+        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")
         {
 
-            // name of this heat source
-            word sourceName = fvSourceCellIndices_.toc()[idxI];
-
-            // sub dictionary with all parameters for this heat source
-            dictionary sourceSubDict = fvSourceSubDict.subDict(sourceName);
-            // NOTE: here power should be in W. We will evenly divide the power by the 
-            // total volume of the source 
-            scalar power = sourceSubDict.getScalar("power");
-
             // loop over all cell indices for this source and assign the source term
             // ----- first loop, calculate the total volume
             scalar totalV = 0.0;
@@ -117,120 +185,57 @@ 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;
-        }
-    }
-    else
-    {
-        FatalErrorIn("calcFvSourceCells") << "source: " << source0 << " not supported!"
-                                          << "Options are: cylinderAnnulusToCell!"
-                                          << abort(FatalError);
-    }
-
-    fvSource.correctBoundaryConditions();
-}
-
-void DAFvSourceHeatSource::calcFvSourceCellIndices(HashTable<labelList>& fvSourceCellIndices)
-{
-    /*
-    Description:
-        Calculate the lists of cell indices that are within the source space
-        NOTE: we support multiple heat sources.
-    
-    Output:
-        fvSourceCellIndices: Hash table that contains the lists of cell indices that 
-        are within the source space. The hash key is the name of the source. We support
-        multiple sources. An example of fvSourceCellIndices can be:
-
-        fvSourceCellIndices=
-        {
-            "source1": {1,2,3,4,5},
-            "source2": {6,7,8,9,10,11,12}
+            Info << "Total heat source for " << sourceName << ": " << sourceTotal << " W." << endl;
         }
-    */
-
-    const dictionary& allOptions = daOption_.getAllOptions();
-
-    dictionary fvSourceSubDict = allOptions.subDict("fvSource");
-
-    if (fvSourceSubDict.toc().size() == 0)
-    {
-        FatalErrorIn("calcSourceCells") << "heatSource is selected as fvSource "
-                                        << " but the options are empty!"
-                                        << abort(FatalError);
-    }
-
-    forAll(fvSourceSubDict.toc(), idxI)
-    {
-        word sourceName = fvSourceSubDict.toc()[idxI];
-
-        fvSourceCellIndices.set(sourceName, {});
-
-        dictionary sourceSubDict = fvSourceSubDict.subDict(sourceName);
-        word sourceType = sourceSubDict.getWord("source");
-        // all available source type are in src/meshTools/sets/cellSources
-        // Example of IO parameters os in applications/utilities/mesh/manipulation/topoSet
-        if (sourceType == "cylinderAnnulusToCell")
+        else if (sourceType == "cylinder")
         {
-            // create a topoSet
-            autoPtr<topoSet> currentSet(
-                topoSet::New(
-                    "cellSet",
-                    mesh_,
-                    "set0",
-                    IOobject::NO_READ));
-            // we need to change the min and max because they need to
-            // be of type point; however, we can't parse point type
-            // in pyDict, we need to change them here.
-
-            scalarList point1;
-            scalarList point2;
-            sourceSubDict.readEntry<scalarList>("p1", point1);
-            sourceSubDict.readEntry<scalarList>("p2", point2);
-
-            point p1;
-            point p2;
-            p1[0] = point1[0];
-            p1[1] = point1[1];
-            p1[2] = point1[2];
-            p2[0] = point2[0];
-            p2[1] = point2[1];
-            p2[2] = point2[2];
-
-            scalar outerRadius = sourceSubDict.getScalar("outerRadius");
-            scalar innerRadius = sourceSubDict.getScalar("innerRadius");
-
-            dictionary tmpDict;
-            tmpDict.set("p1", p1);
-            tmpDict.set("p2", p2);
-            tmpDict.set("innerRadius", innerRadius);
-            tmpDict.set("outerRadius", outerRadius);
-
-            // create the source
-            autoPtr<topoSetSource> sourceSet(
-                topoSetSource::New(sourceType, mesh_, tmpDict));
-
-            // add the sourceSet to topoSet
-            sourceSet().applyToSet(topoSetSource::NEW, currentSet());
-            // get the face index from currentSet, we need to use
-            // this special for loop
-            for (const label i : currentSet())
+            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;
+
+            forAll(mesh_.cells(), cellI)
             {
-                fvSourceCellIndices[sourceName].append(i);
+                // 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);
+
+                if (cellC2AVecRLen <= radius && cellC2AVecALen <= cylinderLen / 2.0)
+                {
+                    fvSource[cellI] += power / totalVolume;
+                }
             }
         }
         else
         {
-            FatalErrorIn("calcFvSourceCells") << "source: " << sourceType << " not supported!"
-                                              << "Options are: cylinderAnnulusToCell!"
-                                              << abort(FatalError);
+            FatalErrorIn("DAFvSourceHeatSource") << "source: " << sourceType << " not supported!"
+                                                 << "Options are: cylinderAnnulusToCell and cylinder!"
+                                                 << abort(FatalError);
         }
     }
 
-    if (daOption_.getOption<label>("debug"))
-    {
-        Info << "fvSourceCellIndices " << fvSourceCellIndices << endl;
-    }
+    fvSource.correctBoundaryConditions();
 }
 
 } // End namespace Foam
diff --git a/src/adjoint/DAFvSource/DAFvSourceHeatSource.H b/src/adjoint/DAFvSource/DAFvSourceHeatSource.H
index 228ad735..236106cd 100755
--- a/src/adjoint/DAFvSource/DAFvSourceHeatSource.H
+++ b/src/adjoint/DAFvSource/DAFvSourceHeatSource.H
@@ -31,8 +31,13 @@ protected:
     /// HashTable that contains lists of cell indices that are within the actuator disk space
     HashTable<labelList> fvSourceCellIndices_;
 
+    HashTable<scalarList> cylinderP1_;
+    HashTable<scalarList> cylinderP2_;
+    HashTable<scalar> cylinderRadius_;
+    HashTable<scalar> power_;
+
     /// calculate DAFvSourceHeatSource::fvSourceCellIndices_
-    void calcFvSourceCellIndices(HashTable<labelList>& fvSourceCellIndices);
+    void calcFvSourceCellIndices(const word sourceName, HashTable<labelList>& fvSourceCellIndices);
 
 public:
     TypeName("heatSource");
diff --git a/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt b/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt
index 0808b305..e96e9016 100755
--- a/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt
+++ b/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt
@@ -1,8 +1,8 @@
 Dictionary Key: HF_INNER
-@value    189999.9999999991 1e-08 1e-10
+@value    181261.4809205142 1e-08 1e-10
 Dictionary Key: HF_OUTER
-@value   -209999.9999999985 1e-08 1e-10
+@value   -218738.5190794858 1e-08 1e-10
 Dictionary Key: TNormSum
-@value    6400011.348643273 1e-08 1e-10
+@value    6400011.350827799 1e-08 1e-10
 Dictionary Key: fail
 @value                    0 1e-08 1e-10
\ No newline at end of file
diff --git a/tests/runTests_DAHeatTransferFoam.py b/tests/runTests_DAHeatTransferFoam.py
index d1d9169a..bc3fb2b8 100755
--- a/tests/runTests_DAHeatTransferFoam.py
+++ b/tests/runTests_DAHeatTransferFoam.py
@@ -53,11 +53,18 @@
     "fvSource": {
         "source1": {
             "type": "heatSource",
-            "source": "cylinderAnnulusToCell",
-            "p1": [0.0, 0.0, 0.0],
-            "p2": [0.5, 0.0, 0.0],
-            "innerRadius": 0.0,
-            "outerRadius": 10.0,
+            "source": "cylinderToCell",
+            "p1": [0.6, 0.055, 0.025],
+            "p2": [1.0, 0.055, 0.025],
+            "radius": 0.01,
+            "power": 1000.0,
+        },
+        "source2": {
+            "type": "heatSource",
+            "source": "cylinder",
+            "p1": [0.0, 0.055, 0.025],
+            "p2": [0.4, 0.055, 0.025],
+            "radius": 0.005,
             "power": 1000.0,
         },
     },