Added an option to snap heatSource center to a cell (#657)

* Fixed an issue for snappedCenter. * Added snapCenterToCell for heat source.
mdolab · Jun 28, 2024 · c195f8b · c195f8b
1 parent a941eaf
commit c195f8b
Show file tree

Hide file tree

Showing 8 changed files with 117 additions and 17 deletions.
diff --git a/src/adjoint/DAFvSource/DAFvSource.C b/src/adjoint/DAFvSource/DAFvSource.C
@@ -199,6 +199,29 @@ void DAFvSource::syncDAOptionToActuatorDVs()
         }
     }
 }
+
+void DAFvSource::findGlobalSnappedCenter(
+    label snappedCenterCellI,
+    vector& center)
+{
+    scalar centerX = 0.0;
+    scalar centerY = 0.0;
+    scalar centerZ = 0.0;
+
+    if (snappedCenterCellI >= 0)
+    {
+        centerX = mesh_.C()[snappedCenterCellI][0];
+        centerY = mesh_.C()[snappedCenterCellI][1];
+        centerZ = mesh_.C()[snappedCenterCellI][2];
+    }
+    reduce(centerX, sumOp<scalar>());
+    reduce(centerY, sumOp<scalar>());
+    reduce(centerZ, sumOp<scalar>());
+
+    center[0] = centerX;
+    center[1] = centerY;
+    center[2] = centerZ;
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

diff --git a/src/adjoint/DAFvSource/DAFvSource.H b/src/adjoint/DAFvSource/DAFvSource.H
@@ -154,6 +154,10 @@ public:
 
     /// virtual function for regIOobject
     bool writeData(Ostream& os) const;
+
+    void findGlobalSnappedCenter(
+        label snappedCenterCellI,
+        vector& center);
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

diff --git a/src/adjoint/DAFvSource/DAFvSourceHeatSource.C b/src/adjoint/DAFvSource/DAFvSourceHeatSource.C
@@ -26,6 +26,12 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
 {
     printInterval_ = daOption.getOption<label>("printInterval");
 
+    // now we need to initialize actuatorDiskDVs_ by synchronizing the values
+    // defined in fvSource from DAOption to actuatorDiskDVs_
+    // NOTE: we need to call this function whenever we change the actuator
+    // design variables during optimization
+    this->syncDAOptionToActuatorDVs();
+
     const dictionary& allOptions = daOption_.getAllOptions();
 
     dictionary fvSourceSubDict = allOptions.subDict("fvSource");
@@ -101,6 +107,28 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
             // eps is a smoothing parameter, it should be the local mesh cell size in meters
             // near the cylinder region
             cylinderEps_.set(sourceName, sourceSubDict.getScalar("eps"));
+
+            snapCenter2Cell_.set(sourceName, sourceSubDict.lookupOrDefault<label>("snapCenter2Cell", 0));
+            if (snapCenter2Cell_[sourceName])
+            {
+                point centerPoint = {actuatorDiskDVs_[sourceName][0], actuatorDiskDVs_[sourceName][1], actuatorDiskDVs_[sourceName][2]};
+                snappedCenterCellI_.set(sourceName, mesh_.findCell(centerPoint));
+                label foundCellI = 0;
+                if (snappedCenterCellI_[sourceName] >= 0)
+                {
+                    foundCellI = 1;
+                }
+                reduce(foundCellI, sumOp<label>());
+                if (foundCellI != 1)
+                {
+                    FatalErrorIn(" ") << "There should be only one cell found globally while "
+                                      << foundCellI << " was returned."
+                                      << " Please adjust center such that it is located completely"
+                                      << " within a cell in the mesh domain. The center should not "
+                                      << " be outside of the mesh domain or on a mesh face "
+                                      << abort(FatalError);
+                }
+            }
         }
         else
         {
@@ -109,12 +137,6 @@ DAFvSourceHeatSource::DAFvSourceHeatSource(
                                                  << abort(FatalError);
         }
     }
-
-    // now we need to initialize actuatorDiskDVs_ by synchronizing the values
-    // defined in fvSource from DAOption to actuatorDiskDVs_
-    // NOTE: we need to call this function whenever we change the actuator
-    // design variables during optimization
-    this->syncDAOptionToActuatorDVs();
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -210,6 +232,11 @@ void DAFvSourceHeatSource::calcFvSource(volScalarField& fvSource)
         {
             vector cylinderCenter =
                 {actuatorDiskDVs_[sourceName][0], actuatorDiskDVs_[sourceName][1], actuatorDiskDVs_[sourceName][2]};
+
+            if (snapCenter2Cell_[sourceName])
+            {
+                this->findGlobalSnappedCenter(snappedCenterCellI_[sourceName], cylinderCenter);
+            }
             vector cylinderDir =
                 {actuatorDiskDVs_[sourceName][3], actuatorDiskDVs_[sourceName][4], actuatorDiskDVs_[sourceName][5]};
             scalar radius = actuatorDiskDVs_[sourceName][6];

diff --git a/src/adjoint/DAFvSource/DAFvSourceHeatSource.H b/src/adjoint/DAFvSource/DAFvSourceHeatSource.H
@@ -40,6 +40,12 @@ protected:
     /// print interval for primal and adjoint
     label printInterval_;
 
+    /// whether to snap the center to a cell in the mesh if yes the center will move with the mesh
+    HashTable<label> snapCenter2Cell_;
+
+    /// the cell index for the center if snapCenter2Cell_ = 1
+    HashTable<label> snappedCenterCellI_;
+
 public:
     TypeName("heatSource");
     // Constructors

diff --git a/src/adjoint/DAObjFunc/DAObjFuncLocation.C b/src/adjoint/DAObjFunc/DAObjFuncLocation.C
@@ -70,9 +70,20 @@ DAObjFuncLocation::DAObjFuncLocation(
     {
         point centerPoint = {center_[0], center_[1], center_[2]};
         snappedCenterCellI_ = mesh_.findCell(centerPoint);
-        if (snappedCenterCellI_ < 0)
+        label foundCellI = 0;
+        if (snappedCenterCellI_ >= 0)
         {
-            FatalErrorIn(" ") << "can not find cell for center " << centerPoint << abort(FatalError);
+            foundCellI = 1;
+        }
+        reduce(foundCellI, sumOp<label>());
+        if (foundCellI != 1)
+        {
+            FatalErrorIn(" ") << "There should be only one cell found globally while "
+                              << foundCellI << " was returned."
+                              << " Please adjust center such that it is located completely"
+                              << " within a cell in the mesh domain. The center should not "
+                              << " be outside of the mesh domain or on a mesh face "
+                              << abort(FatalError);
         }
     }
 
@@ -122,6 +133,29 @@ DAObjFuncLocation::DAObjFuncLocation(
     }
 }
 
+void DAObjFuncLocation::findGlobalSnappedCenter(
+    label snappedCenterCellI,
+    vector& center)
+{
+    scalar centerX = 0.0;
+    scalar centerY = 0.0;
+    scalar centerZ = 0.0;
+
+    if (snappedCenterCellI >= 0)
+    {
+        centerX = mesh_.C()[snappedCenterCellI][0];
+        centerY = mesh_.C()[snappedCenterCellI][1];
+        centerZ = mesh_.C()[snappedCenterCellI][2];
+    }
+    reduce(centerX, sumOp<scalar>());
+    reduce(centerY, sumOp<scalar>());
+    reduce(centerZ, sumOp<scalar>());
+
+    center[0] = centerX;
+    center[1] = centerY;
+    center[2] = centerZ;
+}
+
 /// calculate the value of objective function
 void DAObjFuncLocation::calcObjFunc(
     const labelList& objFuncFaceSources,
@@ -172,7 +206,7 @@ void DAObjFuncLocation::calcObjFunc(
             vector center = center_;
             if (snapCenter2Cell_)
             {
-                center = mesh_.C()[snappedCenterCellI_];
+                this->findGlobalSnappedCenter(snappedCenterCellI_, center);
             }
 
             vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center;
@@ -220,7 +254,7 @@ void DAObjFuncLocation::calcObjFunc(
             vector center = center_;
             if (snapCenter2Cell_)
             {
-                center = mesh_.C()[snappedCenterCellI_];
+                this->findGlobalSnappedCenter(snappedCenterCellI_, center);
             }
 
             vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center;
@@ -261,7 +295,7 @@ void DAObjFuncLocation::calcObjFunc(
             vector center = center_;
             if (snapCenter2Cell_)
             {
-                center = mesh_.C()[snappedCenterCellI_];
+                this->findGlobalSnappedCenter(snappedCenterCellI_, center);
             }
 
             vector faceC = mesh_.Cf().boundaryField()[maxRPatchI_][maxRFaceI_] - center;

diff --git a/src/adjoint/DAObjFunc/DAObjFuncLocation.H b/src/adjoint/DAObjFunc/DAObjFuncLocation.H
@@ -79,6 +79,10 @@ public:
         scalarList& objFuncFaceValues,
         scalarList& objFuncCellValues,
         scalar& objFuncValue);
+
+    void findGlobalSnappedCenter(
+        label snappedCenterCellI,
+        vector& center);
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

diff --git a/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt b/tests/refs/DAFoam_Test_DAHeatTransferFoamRef.txt
@@ -1,14 +1,14 @@
 Dictionary Key: HF_INNER
-@value    179999.9999999994 1e-08 1e-10
+@value     180150.558627527 1e-08 1e-10
 Dictionary Key: HF_OUTER
-@value   -220000.0000000005 1e-08 1e-10
+@value   -220150.5586275289 1e-08 1e-10
 Dictionary Key: IR1KS
 @value    19.89626068314585 1e-08 1e-10
 Dictionary Key: R1
 @value    0.994617581067802 1e-08 1e-10
 Dictionary Key: R1KS
-@value    6.267981796909868 1e-08 1e-10
+@value    6.252992336157548 1e-08 1e-10
 Dictionary Key: TNormSum
-@value     6400011.35124065 1e-08 1e-10
+@value    6400011.351333449 1e-08 1e-10
 Dictionary Key: fail
 @value                    0 1e-08 1e-10
diff --git a/tests/runTests_DAHeatTransferFoam.py b/tests/runTests_DAHeatTransferFoam.py
@@ -59,9 +59,10 @@
                 "patches": ["channel_inner"],
                 "mode": "maxRadiusKS",
                 "axis": [0.0, 1.0, 0.0],
-                "center": [0.5, 0.05, -0.05],
+                "center": [0.505, 0.056, 0.026],
                 "coeffKS": 1.0,
                 "scale": 1.0,
+                "snapCenter2Cell": True,
                 "addToAdjoint": True,
             }
         },
@@ -101,12 +102,13 @@
         "source2": {
             "type": "heatSource",
             "source": "cylinderSmooth",
-            "center": [0.2, 0.055, 0.025],
+            "center": [0.201, 0.056, 0.026],
             "axis": [1.0, 0.0, 0.0],
             "length": 0.4,
             "radius": 0.005,
             "power": 1000.0,
             "eps": 0.001,
+            "snapCenter2Cell": True,
         },
     },
 }