diff --git a/dafoam/mphys/mphys_dafoam.py b/dafoam/mphys/mphys_dafoam.py
index 51ce0be2..f2159379 100644
--- a/dafoam/mphys/mphys_dafoam.py
+++ b/dafoam/mphys/mphys_dafoam.py
@@ -702,7 +702,12 @@ def apply_options(self, optionDict):
     # calculate the residual
     def apply_nonlinear(self, inputs, outputs, residuals):
         DASolver = self.DASolver
-        DASolver.setStates(outputs["%s_states" % self.discipline])
+        # NOTE: we do not pass the states from inputs to the OF layer.
+        # this can cause potential convergence issue because the initial states
+        # in the inputs are set to all ones. So passing this all-ones states
+        # into the OF layer may diverge the primal solver. Here we can always
+        # use the states from the OF layer to compute the residuals.
+        # DASolver.setStates(outputs["%s_states" % self.discipline])
 
         # get flow residuals from DASolver
         residuals["%s_states" % self.discipline] = DASolver.getResiduals()
@@ -1027,6 +1032,9 @@ def solve_linear(self, d_outputs, d_residuals, mode):
                 # in the next line
                 if not self.DASolver.getOption("adjEqnOption")["useNonZeroInitGuess"]:
                     self.psi.set(0)
+                else:
+                    # if useNonZeroInitGuess is True, we will assign the OM's psi to self.psi
+                    self.psi = DASolver.array2Vec(d_residuals["%s_states" % self.discipline].copy())
 
                 if self.DASolver.getOption("adjEqnOption")["dynAdjustTol"]:
                     # if we want to dynamically adjust the tolerance, call this function. This is mostly used
diff --git a/src/adjoint/DAFvSource/DAFvSource.C b/src/adjoint/DAFvSource/DAFvSource.C
index 41901680..a436ebe0 100755
--- a/src/adjoint/DAFvSource/DAFvSource.C
+++ b/src/adjoint/DAFvSource/DAFvSource.C
@@ -100,6 +100,19 @@ void DAFvSource::calcFvSource(volVectorField& fvSource)
                      << abort(FatalError);
 }
 
+void DAFvSource::calcFvSource(volScalarField& fvSource)
+{
+    /*
+    Description:
+        Calculate the fvSource term
+        NOTE: this need to be implemented in the child class, if not,
+        print an error!
+    */
+    FatalErrorIn("") << "calcFvSource not implemented " << endl
+                     << " in the child class for " << modelType_
+                     << abort(FatalError);
+}
+
 bool DAFvSource::writeData(Ostream& os) const
 {
     /*
diff --git a/src/adjoint/DAFvSource/DAFvSource.H b/src/adjoint/DAFvSource/DAFvSource.H
index b5b71d21..e693f125 100755
--- a/src/adjoint/DAFvSource/DAFvSource.H
+++ b/src/adjoint/DAFvSource/DAFvSource.H
@@ -105,6 +105,9 @@ public:
     /// compute the FvSource term
     virtual void calcFvSource(volVectorField& fvSource);
 
+    /// overloaded function for scalar fields
+    virtual void calcFvSource(volScalarField& fvSource);
+
     /// set a new value to the actuator disk design variable
     void setActuatorDVs(
         const word diskName,
diff --git a/src/adjoint/DAFvSource/DAFvSourceHeatSource.C b/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
new file mode 100755
index 00000000..13ec0fa8
--- /dev/null
+++ b/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
@@ -0,0 +1,238 @@
+/*---------------------------------------------------------------------------*\
+
+    DAFoam  : Discrete Adjoint with OpenFOAM
+    Version : v3
+
+\*---------------------------------------------------------------------------*/
+
+#include "DAFvSourceHeatSource.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+defineTypeNameAndDebug(DAFvSourceHeatSource, 0);
+addToRunTimeSelectionTable(DAFvSource, DAFvSourceHeatSource, dictionary);
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+DAFvSourceHeatSource::DAFvSourceHeatSource(
+    const word modelType,
+    const fvMesh& mesh,
+    const DAOption& daOption,
+    const DAModel& daModel,
+    const DAIndex& daIndex)
+    : DAFvSource(modelType, mesh, daOption, daModel, daIndex)
+{
+    this->calcFvSourceCellIndices(fvSourceCellIndices_);
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)
+{
+    /*
+    Description:
+        Calculate the heat source and assign them to fvSource
+
+    Example:
+        An example of the fvSource in pyOptions in pyDAFoam can be
+        defOptions = 
+        {
+            "fvSource"
+            {
+                "source1"
+                {
+                    "type": "heatSource",
+                    "source": "cylinderAnnulusToCell",
+                    "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,
+                    "power": 101.0,  # here we should prescribe the power in W
+                },
+                "source2"
+                {
+                    "type": "heatSource",
+                    "source": "cylinderAnnulusToCell",
+                    "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,
+                    "power": 10.5, # here we should prescribe the power in W
+                }
+            }
+        }
+    */
+
+    forAll(fvSource, idxI)
+    {
+        fvSource[idxI] = 0.0;
+    }
+
+    const dictionary& allOptions = daOption_.getAllOptions();
+
+    dictionary fvSourceSubDict = allOptions.subDict("fvSource");
+
+    word sourceName0 = fvSourceSubDict.toc()[0];
+    word source0 = fvSourceSubDict.subDict(sourceName0).getWord("source");
+
+    if (source0 == "cylinderAnnulusToCell")
+    {
+
+        // loop over all the cell indices for all heat sources
+        forAll(fvSourceCellIndices_.toc(), idxI)
+        {
+
+            // 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;
+            forAll(fvSourceCellIndices_[sourceName], idxJ)
+            {
+                label cellI = fvSourceCellIndices_[sourceName][idxJ];
+                scalar cellV = mesh_.V()[cellI];
+                totalV += cellV;
+            }
+            reduce(totalV, sumOp<scalar>());
+
+            // ------- second loop, assign power
+            scalar sourceTotal = 0.0;
+            forAll(fvSourceCellIndices_[sourceName], idxJ)
+            {
+                // cell index
+                label cellI = fvSourceCellIndices_[sourceName][idxJ];
+
+                fvSource[cellI] = power / totalV;
+                sourceTotal += fvSource[cellI] * mesh_.V()[cellI];
+            }
+
+            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}
+        }
+    */
+
+    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")
+        {
+            // 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())
+            {
+                fvSourceCellIndices[sourceName].append(i);
+            }
+        }
+        else
+        {
+            FatalErrorIn("calcFvSourceCells") << "source: " << sourceType << " not supported!"
+                                              << "Options are: cylinderAnnulusToCell!"
+                                              << abort(FatalError);
+        }
+    }
+
+    if (daOption_.getOption<label>("debug"))
+    {
+        Info << "fvSourceCellIndices " << fvSourceCellIndices << endl;
+    }
+}
+
+} // End namespace Foam
+
+// ************************************************************************* //
diff --git a/src/adjoint/DAFvSource/DAFvSourceHeatSource.H b/src/adjoint/DAFvSource/DAFvSourceHeatSource.H
new file mode 100755
index 00000000..228ad735
--- /dev/null
+++ b/src/adjoint/DAFvSource/DAFvSourceHeatSource.H
@@ -0,0 +1,66 @@
+/*---------------------------------------------------------------------------*\
+
+    DAFoam  : Discrete Adjoint with OpenFOAM
+    Version : v3
+
+    Description:
+        Child class for HeatSource disk source
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DAFvSourceHeatSource_H
+#define DAFvSourceHeatSource_H
+
+#include "DAFvSource.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+      Class DAFvSourceHeatSource Declaration
+\*---------------------------------------------------------------------------*/
+
+class DAFvSourceHeatSource
+    : public DAFvSource
+{
+
+protected:
+    /// HashTable that contains lists of cell indices that are within the actuator disk space
+    HashTable<labelList> fvSourceCellIndices_;
+
+    /// calculate DAFvSourceHeatSource::fvSourceCellIndices_
+    void calcFvSourceCellIndices(HashTable<labelList>& fvSourceCellIndices);
+
+public:
+    TypeName("heatSource");
+    // Constructors
+
+    //- Construct from components
+    DAFvSourceHeatSource(
+        const word modelType,
+        const fvMesh& mesh,
+        const DAOption& daOption,
+        const DAModel& daModel,
+        const DAIndex& daIndex);
+
+    //- Destructor
+    virtual ~DAFvSourceHeatSource()
+    {
+    }
+
+    /// compute the FvSource term
+    virtual void calcFvSource(volScalarField& fvSource);
+};
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.C b/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.C
index bd6595d6..f64a2e2d 100755
--- a/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.C
+++ b/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.C
@@ -22,7 +22,8 @@ DAHeatTransferFoam::DAHeatTransferFoam(
     : DASolver(argsAll, pyOptions),
       TPtr_(nullptr),
       fvSourcePtr_(nullptr),
-      kPtr_(nullptr)
+      kPtr_(nullptr),
+      daFvSourcePtr_(nullptr)
 {
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -45,6 +46,22 @@ void DAHeatTransferFoam::initSolver()
     daLinearEqnPtr_.reset(new DALinearEqn(mesh, daOptionPtr_()));
 
     this->setDAObjFuncList();
+
+    // initialize fvSource and compute the source term
+    const dictionary& allOptions = daOptionPtr_->getAllOptions();
+    if (allOptions.subDict("fvSource").toc().size() != 0)
+    {
+        hasFvSource_ = 1;
+        Info << "Initializing DASource" << endl;
+        word sourceName = allOptions.subDict("fvSource").toc()[0];
+        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);
+    }
 }
 
 label DAHeatTransferFoam::solvePrimal(
@@ -82,7 +99,7 @@ label DAHeatTransferFoam::solvePrimal(
         return 1;
     }
 
-    label printInterval = daOptionPtr_->getOption<label>("printIntervalUnsteady");
+    label printInterval = daOptionPtr_->getOption<label>("printInterval");
     label printToScreen = 0;
     // main loop
     while (this->loop(runTime)) // using simple.loop() will have seg fault in parallel
diff --git a/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.H b/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.H
index 200d4406..33cbd81d 100755
--- a/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.H
+++ b/src/adjoint/DASolver/DAHeatTransferFoam/DAHeatTransferFoam.H
@@ -12,6 +12,7 @@
 #define DAHeatTransferFoam_H
 
 #include "DASolver.H"
+#include "DAFvSource.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -37,6 +38,12 @@ protected:
     /// h pointer
     autoPtr<dimensionedScalar> kPtr_;
 
+    /// DASource pointer
+    autoPtr<DAFvSource> daFvSourcePtr_;
+
+    /// whether to have fvSource term
+    label hasFvSource_ = 0;
+
 public:
     TypeName("DAHeatTransferFoam");
     // Constructors
diff --git a/src/adjoint/Make/files_Solid b/src/adjoint/Make/files_Solid
index eaccdd0e..e4d9da1e 100755
--- a/src/adjoint/Make/files_Solid
+++ b/src/adjoint/Make/files_Solid
@@ -29,6 +29,7 @@ DAObjFunc/DAObjFuncWallHeatFlux.C
 DAObjFunc/DAObjFuncPatchMean.C
 
 DAFvSource/DAFvSource.C
+DAFvSource/DAFvSourceHeatSource.C
 
 DAResidual/DAResidual.C
 DAResidual/DAResidualDummy.C
diff --git a/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt b/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt
index 3310313a..0808b305 100755
--- a/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt
+++ b/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt
@@ -1,8 +1,8 @@
 Dictionary Key: HF_INNER
-@value    199999.9999999992 1e-08 1e-10
+@value    189999.9999999991 1e-08 1e-10
 Dictionary Key: HF_OUTER
-@value   -199999.9999999997 1e-08 1e-10
+@value   -209999.9999999985 1e-08 1e-10
 Dictionary Key: TNormSum
-@value    6400011.343935275 1e-08 1e-10
+@value    6400011.348643273 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 254ad8e4..d1d9169a 100755
--- a/tests/runTests_DAHeatTransferFoam.py
+++ b/tests/runTests_DAHeatTransferFoam.py
@@ -50,6 +50,17 @@
             },
         }
     },
+    "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,
+            "power": 1000.0,
+        },
+    },
 }
 DASolver = PYDAFOAM(options=aeroOptions, comm=MPI.COMM_WORLD)