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] 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,