diff --git a/src/adjoint/DAFvSource/DAFvSourceHeatSource.C b/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
index f86b5fb0..a3a7261d 100755
--- a/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
+++ b/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
@@ -112,11 +112,16 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
             if (snapCenter2Cell_[sourceName])
             {
                 point centerPoint = {actuatorDiskDVs_[sourceName][0], actuatorDiskDVs_[sourceName][1], actuatorDiskDVs_[sourceName][2]};
-                snappedCenterCellI_.set(sourceName, mesh_.findCell(centerPoint));
+
+                // NOTE: we need to call a self-defined findCell func to make it work correctly in ADR
+                label myCellI = DAUtility::myFindCell(mesh_, centerPoint);
+
+                snappedCenterCellI_.set(sourceName, myCellI);
                 label foundCellI = 0;
                 if (snappedCenterCellI_[sourceName] >= 0)
                 {
                     foundCellI = 1;
+                    //Pout << "snap source " << sourceName << " to center " << mesh_.C()[snappedCenterCellI_[sourceName]] << endl;
                 }
                 reduce(foundCellI, sumOp<label>());
                 if (foundCellI != 1)
@@ -128,6 +133,10 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
                                       << " be outside of the mesh domain or on a mesh face "
                                       << abort(FatalError);
                 }
+
+                vector snappedCenter = vector::zero;
+                this->findGlobalSnappedCenter(snappedCenterCellI_[sourceName], snappedCenter);
+                Info << "heat source " << sourceName << " snap to center " << snappedCenter << endl;
             }
         }
         else
@@ -232,7 +241,7 @@ void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)
         {
             vector cylinderCenter =
                 {actuatorDiskDVs_[sourceName][0], actuatorDiskDVs_[sourceName][1], actuatorDiskDVs_[sourceName][2]};
-            
+
             if (snapCenter2Cell_[sourceName])
             {
                 this->findGlobalSnappedCenter(snappedCenterCellI_[sourceName], cylinderCenter);
diff --git a/src/adjoint/DAObjFunc/DAObjFuncLocation.C b/src/adjoint/DAObjFunc/DAObjFuncLocation.C
index 58bbb5fa..f10b0ba7 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncLocation.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncLocation.C
@@ -69,7 +69,10 @@ DAObjFuncLocation::DAObjFuncLocation(
     if (snapCenter2Cell_)
     {
         point centerPoint = {center_[0], center_[1], center_[2]};
-        snappedCenterCellI_ = mesh_.findCell(centerPoint);
+
+        // NOTE: we need to call a self-defined findCell func to make it work correctly in ADR
+        snappedCenterCellI_ = DAUtility::myFindCell(mesh_, centerPoint);
+        
         label foundCellI = 0;
         if (snappedCenterCellI_ >= 0)
         {
@@ -85,6 +88,9 @@ DAObjFuncLocation::DAObjFuncLocation(
                               << " be outside of the mesh domain or on a mesh face "
                               << abort(FatalError);
         }
+        vector snappedCenter = vector::zero;
+        this->findGlobalSnappedCenter(snappedCenterCellI_, snappedCenter);
+        Info << "snap to center " << snappedCenter << endl;
     }
 
     if (mode_ == "maxRadius")
@@ -196,6 +202,12 @@ void DAObjFuncLocation::calcObjFunc(
         // calculate Location
         scalar objValTmp = 0.0;
 
+        vector center = center_;
+        if (snapCenter2Cell_)
+        {
+            this->findGlobalSnappedCenter(snappedCenterCellI_, center);
+        }
+
         forAll(objFuncFaceSources, idxI)
         {
             const label& objFuncFaceI = objFuncFaceSources[idxI];
@@ -203,12 +215,6 @@ void DAObjFuncLocation::calcObjFunc(
             const label patchI = daIndex_.bFacePatchI[bFaceI];
             const label faceI = daIndex_.bFaceFaceI[bFaceI];
 
-            vector center = center_;
-            if (snapCenter2Cell_)
-            {
-                this->findGlobalSnappedCenter(snappedCenterCellI_, center);
-            }
-
             vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center;
 
             tensor faceCTensor(tensor::zero);
@@ -244,6 +250,12 @@ void DAObjFuncLocation::calcObjFunc(
         // calculate Location
         scalar objValTmp = 0.0;
 
+        vector center = center_;
+        if (snapCenter2Cell_)
+        {
+            this->findGlobalSnappedCenter(snappedCenterCellI_, center);
+        }
+
         forAll(objFuncFaceSources, idxI)
         {
             const label& objFuncFaceI = objFuncFaceSources[idxI];
@@ -251,12 +263,6 @@ void DAObjFuncLocation::calcObjFunc(
             const label patchI = daIndex_.bFacePatchI[bFaceI];
             const label faceI = daIndex_.bFaceFaceI[bFaceI];
 
-            vector center = center_;
-            if (snapCenter2Cell_)
-            {
-                this->findGlobalSnappedCenter(snappedCenterCellI_, center);
-            }
-
             vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center;
 
             tensor faceCTensor(tensor::zero);
@@ -289,14 +295,15 @@ void DAObjFuncLocation::calcObjFunc(
     else if (mode_ == "maxRadius")
     {
         scalar radius = 0.0;
-        if (maxRPatchI_ >= 0 && maxRFaceI_ >= 0)
+
+        vector center = center_;
+        if (snapCenter2Cell_)
         {
+            this->findGlobalSnappedCenter(snappedCenterCellI_, center);
+        }
 
-            vector center = center_;
-            if (snapCenter2Cell_)
-            {
-                this->findGlobalSnappedCenter(snappedCenterCellI_, center);
-            }
+        if (maxRPatchI_ >= 0 && maxRFaceI_ >= 0)
+        {
 
             vector faceC = mesh_.Cf().boundaryField()[maxRPatchI_][maxRFaceI_] - center;
 
diff --git a/src/adjoint/DAObjFunc/DAObjFuncMass.C b/src/adjoint/DAObjFunc/DAObjFuncMass.C
index 053a9941..f40e804e 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncMass.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncMass.C
@@ -44,6 +44,7 @@ DAObjFuncMass::DAObjFuncMass(
 
     objFuncDict_.readEntry<scalar>("scale", scale_);
 
+    rho_ = objFuncDict_.lookupOrDefault<scalar>("rho", -1.0);
 }
 
 /// calculate the value of objective function
@@ -82,15 +83,30 @@ void DAObjFuncMass::calcObjFunc(
     objFuncValue = 0.0;
 
     const objectRegistry& db = mesh_.thisDb();
-    const volScalarField& rho = db.lookupObject<volScalarField>("solid:rho");
 
-    // calculate mass
-    forAll(objFuncCellSources, idxI)
+    if (rho_ < 0.0)
     {
-        const label& cellI = objFuncCellSources[idxI];
-        scalar volume = mesh_.V()[cellI];
-        objFuncCellValues[idxI] = scale_ * volume * rho[cellI];
-        objFuncValue += objFuncCellValues[idxI];
+        const volScalarField& rho = db.lookupObject<volScalarField>("solid:rho");
+
+        // calculate mass
+        forAll(objFuncCellSources, idxI)
+        {
+            const label& cellI = objFuncCellSources[idxI];
+            scalar volume = mesh_.V()[cellI];
+            objFuncCellValues[idxI] = scale_ * volume * rho[cellI];
+            objFuncValue += objFuncCellValues[idxI];
+        }
+    }
+    else
+    {
+        // calculate mass
+        forAll(objFuncCellSources, idxI)
+        {
+            const label& cellI = objFuncCellSources[idxI];
+            scalar volume = mesh_.V()[cellI];
+            objFuncCellValues[idxI] = scale_ * volume * rho_;
+            objFuncValue += objFuncCellValues[idxI];
+        }
     }
 
     // need to reduce the sum of force across all processors
diff --git a/src/adjoint/DAObjFunc/DAObjFuncMass.H b/src/adjoint/DAObjFunc/DAObjFuncMass.H
index 01d0389e..ccdf29e6 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncMass.H
+++ b/src/adjoint/DAObjFunc/DAObjFuncMass.H
@@ -31,6 +31,9 @@ public:
     TypeName("mass");
     // Constructors
 
+    /// user-defined density
+    scalar rho_ = -1.0;
+
     //- Construct from components
     DAObjFuncMass(
         const fvMesh& mesh,
diff --git a/src/adjoint/DAObjFunc/DAObjFuncVariance.C b/src/adjoint/DAObjFunc/DAObjFuncVariance.C
index adf5633c..a32e40ee 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncVariance.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncVariance.C
@@ -207,7 +207,7 @@ DAObjFuncVariance::DAObjFuncVariance(
                     forAll(probePointCoords_, idxI)
                     {
                         point pointCoord = {probePointCoords_[idxI][0], probePointCoords_[idxI][1], probePointCoords_[idxI][2]};
-                        label cellI = mesh_.findCell(pointCoord);
+                        label cellI = DAUtility::myFindCell(mesh_, pointCoord);
                         if (cellI >= 0)
                         {
                             probeCellIndex_.append(cellI);
@@ -286,7 +286,7 @@ DAObjFuncVariance::DAObjFuncVariance(
                     forAll(probePointCoords_, idxI)
                     {
                         point pointCoord = {probePointCoords_[idxI][0], probePointCoords_[idxI][1], probePointCoords_[idxI][2]};
-                        label cellI = mesh_.findCell(pointCoord);
+                        label cellI = DAUtility::myFindCell(mesh_, pointCoord);
                         if (cellI >= 0)
                         {
                             probeCellIndex_.append(cellI);
diff --git a/src/adjoint/DAResidual/DAResidualHeatTransferFoam.H b/src/adjoint/DAResidual/DAResidualHeatTransferFoam.H
index 4cea3f32..5adca833 100755
--- a/src/adjoint/DAResidual/DAResidualHeatTransferFoam.H
+++ b/src/adjoint/DAResidual/DAResidualHeatTransferFoam.H
@@ -37,7 +37,7 @@ protected:
 
     autoPtr<dimensionedScalar> kPtr_;
 
-    label hasFvSource_;
+    label hasFvSource_ = 0;
 
 public:
     TypeName("DAHeatTransferFoam");
diff --git a/src/adjoint/DAUtility/DAUtility.C b/src/adjoint/DAUtility/DAUtility.C
index 423740e6..e34737fc 100755
--- a/src/adjoint/DAUtility/DAUtility.C
+++ b/src/adjoint/DAUtility/DAUtility.C
@@ -770,7 +770,7 @@ void DAUtility::primalResidualControl(
     // calculate the initial residual mag and set it to primalResidualNorms_
 
     // for vectors, we need to use the median value for the residual
-    // this is because we often need to run 2D simulations with symmetry 
+    // this is because we often need to run 2D simulations with symmetry
     // BC, so one component of the residual vector, which is related to the symmetry BC,
     // may be high while the other two components' residuals are low.
     // In this case, we can use the median value for the residual vector, which better
@@ -792,6 +792,20 @@ void DAUtility::primalResidualControl(
     }
 }
 
+label DAUtility::myFindCell(
+    const primitiveMesh& mesh,
+    const point& point)
+{
+    /*
+    A self-defined findCell function. We will need to cast fvMesh to primitiveMesh
+    and then call primitiveMesh's findCell. For some reasons, the fvMesh's findCell
+    did not work correctly in ADR mode...
+    */
+
+    label cellI = mesh.findCell(point);
+    return cellI;
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/adjoint/DAUtility/DAUtility.H b/src/adjoint/DAUtility/DAUtility.H
index 0ea43c5e..97b4b65e 100755
--- a/src/adjoint/DAUtility/DAUtility.H
+++ b/src/adjoint/DAUtility/DAUtility.H
@@ -144,12 +144,16 @@ public:
         const SolverPerformance<scalar>& solverP,
         const label printToScreen,
         const word varName);
-    
+
     /// control when to print the residual and also compute the maxInitRes
     static void primalResidualControl(
         const SolverPerformance<vector>& solverP,
         const label printToScreen,
         const word varName);
+
+    static label myFindCell(
+        const primitiveMesh& mesh,
+        const point& point);
 };
 
 template<class classType>