Added inverseRadius for localtion obj. (mdolab#654)

friedenhe · Jun 24, 2024 · aa2d351 · aa2d351
1 parent 7b3f157
commit aa2d351
Show file tree

Hide file tree

Showing 2 changed files with 43 additions and 3 deletions.
diff --git a/src/adjoint/DAObjFunc/DAObjFuncLocation.C b/src/adjoint/DAObjFunc/DAObjFuncLocation.C
@@ -109,7 +109,6 @@ DAObjFuncLocation::DAObjFuncLocation(
             maxRFaceI_ = maxRFaceI;
         }
     }
-
 }
 
 /// calculate the value of objective function
@@ -185,8 +184,48 @@ void DAObjFuncLocation::calcObjFunc(
         // need to reduce the sum of force across all processors
         reduce(objValTmp, sumOp<scalar>());
 
-        // expSumKS stores sum[exp(coeffKS*x_i)], it will be used to scale dFdW
-        expSumKS = objValTmp;
+        objFuncValue = log(objValTmp) / coeffKS_;
+    }
+    else if (mode_ == "maxInverseRadiusKS")
+    {
+        // this is essentially minimal radius using KS
+
+        // calculate Location
+        scalar objValTmp = 0.0;
+
+        forAll(objFuncFaceSources, idxI)
+        {
+            const label& objFuncFaceI = objFuncFaceSources[idxI];
+            label bFaceI = objFuncFaceI - daIndex_.nLocalInternalFaces;
+            const label patchI = daIndex_.bFacePatchI[bFaceI];
+            const label faceI = daIndex_.bFaceFaceI[bFaceI];
+
+            vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center_;
+
+            tensor faceCTensor(tensor::zero);
+            faceCTensor.xx() = faceC.x();
+            faceCTensor.yy() = faceC.y();
+            faceCTensor.zz() = faceC.z();
+
+            vector faceCAxial = faceCTensor & axis_;
+            vector faceCRadial = faceC - faceCAxial;
+
+            scalar radius = mag(faceCRadial);
+            scalar iRadius = 1.0 / (radius + 1e-12);
+
+            objFuncFaceValues[idxI] = exp(coeffKS_ * iRadius);
+
+            objValTmp += objFuncFaceValues[idxI];
+
+            if (objValTmp > 1e200)
+            {
+                FatalErrorIn(" ") << "KS function summation term too large! "
+                                  << "Reduce coeffKS! " << abort(FatalError);
+            }
+        }
+
+        // need to reduce the sum of force across all processors
+        reduce(objValTmp, sumOp<scalar>());
 
         objFuncValue = log(objValTmp) / coeffKS_;
     }

diff --git a/src/adjoint/Make/files_Solid b/src/adjoint/Make/files_Solid
@@ -27,6 +27,7 @@ DAObjFunc/DAObjFuncVariableVolSum.C
 DAObjFunc/DAObjFuncMeshQualityKS.C
 DAObjFunc/DAObjFuncWallHeatFlux.C
 DAObjFunc/DAObjFuncPatchMean.C
+DAObjFunc/DAObjFuncLocation.C
 
 DAFvSource/DAFvSource.C
 DAFvSource/DAFvSourceHeatSource.C