From 19f212d0e7e5ea678dce1ed707ef11010c6833dc Mon Sep 17 00:00:00 2001
From: Ping He <phe@iastate.edu>
Date: Sat, 29 Jun 2024 17:28:26 -0500
Subject: [PATCH 1/4] Fixed bugs for snapCenter2Cell and hasFvSource.

---
 src/adjoint/DAFvSource/DAFvSourceHeatSource.C    |  6 +++++-
 src/adjoint/DAObjFunc/DAObjFuncFieldInversion.C  |  2 +-
 src/adjoint/DAObjFunc/DAObjFuncLocation.C        |  5 ++++-
 src/adjoint/DAObjFunc/DAObjFuncVariance.C        |  4 ++--
 .../DAResidual/DAResidualHeatTransferFoam.H      |  2 +-
 src/adjoint/DAUtility/DAUtility.C                | 16 +++++++++++++++-
 src/adjoint/DAUtility/DAUtility.H                |  6 +++++-
 7 files changed, 33 insertions(+), 8 deletions(-)

diff --git a/src/adjoint/DAFvSource/DAFvSourceHeatSource.C b/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
index f86b5fb0..975c1ed3 100755
--- a/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
+++ b/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
@@ -112,7 +112,11 @@ 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)
                 {
diff --git a/src/adjoint/DAObjFunc/DAObjFuncFieldInversion.C b/src/adjoint/DAObjFunc/DAObjFuncFieldInversion.C
index b8608271..4fc8ce80 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncFieldInversion.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncFieldInversion.C
@@ -250,7 +250,7 @@ void DAObjFuncFieldInversion::calcObjFunc(
             pRefCoords_[2] = pRefCoords[2];
 
             pRef_ = 0.0;
-            label cellID = mesh_.findCell(pRefCoords_);
+            label cellID = DAUtility::myFindCell(pRefCoords_);
             // only assign pRef if the required cell is found in processor
             if (cellID != -1)
             {
diff --git a/src/adjoint/DAObjFunc/DAObjFuncLocation.C b/src/adjoint/DAObjFunc/DAObjFuncLocation.C
index 58bbb5fa..b557872e 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)
         {
diff --git a/src/adjoint/DAObjFunc/DAObjFuncVariance.C b/src/adjoint/DAObjFunc/DAObjFuncVariance.C
index adf5633c..5ac06d6c 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(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(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>

From 2236a1069af1d404cca65c765c645d643d8fc39e Mon Sep 17 00:00:00 2001
From: Ping He <phe@iastate.edu>
Date: Sat, 29 Jun 2024 18:13:15 -0500
Subject: [PATCH 2/4] Fixed a bug.

---
 src/adjoint/DAObjFunc/DAObjFuncFieldInversion.C | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/adjoint/DAObjFunc/DAObjFuncFieldInversion.C b/src/adjoint/DAObjFunc/DAObjFuncFieldInversion.C
index 4fc8ce80..b8608271 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncFieldInversion.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncFieldInversion.C
@@ -250,7 +250,7 @@ void DAObjFuncFieldInversion::calcObjFunc(
             pRefCoords_[2] = pRefCoords[2];
 
             pRef_ = 0.0;
-            label cellID = DAUtility::myFindCell(pRefCoords_);
+            label cellID = mesh_.findCell(pRefCoords_);
             // only assign pRef if the required cell is found in processor
             if (cellID != -1)
             {

From 465a738d8cc54dd128e29479d060e152f24c35cc Mon Sep 17 00:00:00 2001
From: Ping He <phe@iastate.edu>
Date: Sat, 29 Jun 2024 18:36:46 -0500
Subject: [PATCH 3/4] More bug fix.

---
 src/adjoint/DAObjFunc/DAObjFuncVariance.C | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/adjoint/DAObjFunc/DAObjFuncVariance.C b/src/adjoint/DAObjFunc/DAObjFuncVariance.C
index 5ac06d6c..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 = DAUtility::myFindCell(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 = DAUtility::myFindCell(pointCoord);
+                        label cellI = DAUtility::myFindCell(mesh_, pointCoord);
                         if (cellI >= 0)
                         {
                             probeCellIndex_.append(cellI);

From ac5ff75699f65be5a5d00dddd23ce12cd31affa5 Mon Sep 17 00:00:00 2001
From: Ping He <phe@iastate.edu>
Date: Sun, 30 Jun 2024 08:48:08 -0500
Subject: [PATCH 4/4] Fixed bugs in snapCenter.

---
 src/adjoint/DAFvSource/DAFvSourceHeatSource.C |  7 +++-
 src/adjoint/DAObjFunc/DAObjFuncLocation.C     | 40 ++++++++++---------
 src/adjoint/DAObjFunc/DAObjFuncMass.C         | 30 ++++++++++----
 src/adjoint/DAObjFunc/DAObjFuncMass.H         |  3 ++
 4 files changed, 54 insertions(+), 26 deletions(-)

diff --git a/src/adjoint/DAFvSource/DAFvSourceHeatSource.C b/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
index 975c1ed3..a3a7261d 100755
--- a/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
+++ b/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
@@ -121,6 +121,7 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
                 if (snappedCenterCellI_[sourceName] >= 0)
                 {
                     foundCellI = 1;
+                    //Pout << "snap source " << sourceName << " to center " << mesh_.C()[snappedCenterCellI_[sourceName]] << endl;
                 }
                 reduce(foundCellI, sumOp<label>());
                 if (foundCellI != 1)
@@ -132,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
@@ -236,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 b557872e..f10b0ba7 100755
--- a/src/adjoint/DAObjFunc/DAObjFuncLocation.C
+++ b/src/adjoint/DAObjFunc/DAObjFuncLocation.C
@@ -88,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")
@@ -199,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];
@@ -206,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);
@@ -247,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];
@@ -254,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);
@@ -292,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,