diff --git a/src/adjoint/DAObjFunc/DAObjFuncLocation.C b/src/adjoint/DAObjFunc/DAObjFuncLocation.C index 828adc43..6152ed57 100755 --- a/src/adjoint/DAObjFunc/DAObjFuncLocation.C +++ b/src/adjoint/DAObjFunc/DAObjFuncLocation.C @@ -65,7 +65,50 @@ DAObjFuncLocation::DAObjFuncLocation( center_ = {centerRead[0], centerRead[1], centerRead[2]}; } - + if (mode_ == "maxRadius") + { + // we need to identify the patchI and faceI that has maxR + // the proc that does not own the maxR face will have negative patchI and faceI + // we assume the patchI and faceI do not change during the optimization + // otherwise, we should use maxRadiusKS instead + scalar maxR = -100000; + label maxRPatchI, maxRFaceI; + 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); + + if (radius > maxR) + { + maxR = radius; + maxRPatchI = patchI; + maxRFaceI = faceI; + } + } + + scalar maxRGlobal = maxR; + reduce(maxRGlobal, maxOp()); + + if (fabs(maxRGlobal - maxR) < 1e-8) + { + maxRPatchI_ = maxRPatchI; + maxRFaceI_ = maxRFaceI; + } + } } /// calculate the value of objective function @@ -103,19 +146,20 @@ void DAObjFuncLocation::calcObjFunc( // initialize objFunValue objFuncValue = 0.0; - // calculate Location - scalar objValTmp = 0.0; - forAll(objFuncFaceSources, idxI) + if (mode_ == "maxRadiusKS") { - const label& objFuncFaceI = objFuncFaceSources[idxI]; - label bFaceI = objFuncFaceI - daIndex_.nLocalInternalFaces; - const label patchI = daIndex_.bFacePatchI[bFaceI]; - const label faceI = daIndex_.bFaceFaceI[bFaceI]; + // calculate Location + scalar objValTmp = 0.0; - vector faceC = mesh_.Cf().boundaryField()[patchI][faceI] - center_; - - if (mode_ == "maxRadius") + 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(); @@ -136,21 +180,43 @@ void DAObjFuncLocation::calcObjFunc( << "Reduce coeffKS! " << abort(FatalError); } } - else - { - FatalErrorIn("DAObjFuncLocation") << "mode: " << mode_ << " not supported!" - << "Options are: maxRadius" - << abort(FatalError); - } + + // need to reduce the sum of force across all processors + reduce(objValTmp, sumOp()); + + // expSumKS stores sum[exp(coeffKS*x_i)], it will be used to scale dFdW + expSumKS = objValTmp; + + objFuncValue = log(objValTmp) / coeffKS_; } + else if (mode_ == "maxRadius") + { + scalar radius = 0.0; + if (maxRPatchI_ >= 0 && maxRFaceI_ >= 0) + { + vector faceC = mesh_.Cf().boundaryField()[maxRPatchI_][maxRFaceI_] - center_; + + tensor faceCTensor(tensor::zero); + faceCTensor.xx() = faceC.x(); + faceCTensor.yy() = faceC.y(); + faceCTensor.zz() = faceC.z(); - // need to reduce the sum of force across all processors - reduce(objValTmp, sumOp()); + vector faceCAxial = faceCTensor & axis_; + vector faceCRadial = faceC - faceCAxial; - // expSumKS stores sum[exp(coeffKS*x_i)], it will be used to scale dFdW - expSumKS = objValTmp; + radius = mag(faceCRadial); + } + + reduce(radius, sumOp()); - objFuncValue = log(objValTmp) / coeffKS_; + objFuncValue = radius; + } + else + { + FatalErrorIn("DAObjFuncLocation") << "mode: " << mode_ << " not supported!" + << "Options are: maxRadius" + << abort(FatalError); + } return; } diff --git a/src/adjoint/DAObjFunc/DAObjFuncLocation.H b/src/adjoint/DAObjFunc/DAObjFuncLocation.H index e2028555..23026e7b 100755 --- a/src/adjoint/DAObjFunc/DAObjFuncLocation.H +++ b/src/adjoint/DAObjFunc/DAObjFuncLocation.H @@ -37,8 +37,15 @@ protected: /// axis for radius computation vector axis_ = {1, 0, 0}; + /// center for radius computation vector center_ = {0, 0, 0}; + /// maxRadius patchI + label maxRPatchI_ = -1; + + /// maxRadius patchI + label maxRFaceI_ = -1; + public: TypeName("location"); // Constructors diff --git a/tests/refs/DAFoam_Test_DASimpleFoamADRef.txt b/tests/refs/DAFoam_Test_DASimpleFoamADRef.txt index 32b5312c..9386c0b7 100755 --- a/tests/refs/DAFoam_Test_DASimpleFoamADRef.txt +++ b/tests/refs/DAFoam_Test_DASimpleFoamADRef.txt @@ -5,6 +5,8 @@ Dictionary Key: CL Dictionary Key: COP @value 0.2298218575123644 1e-08 1e-10 Dictionary Key: R1 +@value 0.9988274383772859 1e-08 1e-10 +Dictionary Key: R1KS @value 1.022285515682389 1e-08 1e-10 Dictionary Key: fail @value 0 1e-08 1e-10 @@ -17,6 +19,8 @@ Dictionary Key: shapey @value -0.01063654962455561 1e-05 1e-07 @value 0.005827183649506049 1e-05 1e-07 @value 0.01999380242276746 1e-05 1e-07 +Dictionary Key: twist +@value 0.002313367328235666 1e-05 1e-07 Dictionary Key: CL Dictionary Key: alpha @value 0.08995513866374753 1e-05 1e-07 @@ -24,6 +28,8 @@ Dictionary Key: shapey @value 0.7008636073013889 1e-05 1e-07 @value 0.6091514298213649 1e-05 1e-07 @value 1.418731716621286 1e-05 1e-07 +Dictionary Key: twist +@value 0.09008667655578639 1e-05 1e-07 Dictionary Key: COP Dictionary Key: alpha @value -0.002070288134429892 1e-05 1e-07 @@ -31,13 +37,26 @@ Dictionary Key: shapey @value 0.02508770027211096 1e-05 1e-07 @value 0.397331785797213 1e-05 1e-07 @value 1.137565574218442 1e-05 1e-07 +Dictionary Key: twist +@value -0.001378560585184985 1e-05 1e-07 Dictionary Key: R1 Dictionary Key: alpha @value 0 1e-05 1e-07 Dictionary Key: shapey +@value -2.09670026639443e-10 1e-05 1e-07 +@value -8.466467024008695e-08 1e-05 1e-07 +@value -1.136770130241234e-05 1e-05 1e-07 +Dictionary Key: twist +@value -4.573679551391674e-06 1e-05 1e-07 +Dictionary Key: R1KS +Dictionary Key: alpha +@value 0 1e-05 1e-07 +Dictionary Key: shapey @value 3.445986662700833e-10 1e-05 1e-07 @value 3.454954650172768e-08 1e-05 1e-07 @value 1.51360700044254e-06 1e-05 1e-07 +Dictionary Key: twist +@value -1.990527426047606e-20 1e-05 1e-07 Dictionary Key: dForcedW @value 1.968811128068832 1e-05 1e-07 Dictionary Key: dForcedxV diff --git a/tests/runTests_DASimpleFoamAD.py b/tests/runTests_DASimpleFoamAD.py index 60a3fc5f..aace24ed 100755 --- a/tests/runTests_DASimpleFoamAD.py +++ b/tests/runTests_DASimpleFoamAD.py @@ -94,11 +94,22 @@ "patches": ["wing"], "mode": "maxRadius", "axis": [0.0, 0.0, 1.0], + "scale": 1.0, + "addToAdjoint": True, + } + }, + "R1KS": { + "part1": { + "type": "location", + "source": "patchToFace", + "patches": ["wing"], + "mode": "maxRadiusKS", + "axis": [0.0, 0.0, 1.0], "coeffKS": 100.0, "scale": 1.0, "addToAdjoint": True, } - } + }, }, "normalizeStates": {"U": U0, "p": U0 * U0 / 2.0, "k": k0, "omega": omega0, "phi": 1.0}, "adjPartDerivFDStep": {"State": 1e-6, "FFD": 1e-3, "ACTD": 1.0e-3}, @@ -107,6 +118,7 @@ "designVar": { "shapey": {"designVarType": "FFD"}, "alpha": {"designVarType": "AOA", "patches": ["inout"], "flowAxis": "x", "normalAxis": "y"}, + "twist": {"designVarType": "FFD"} }, } @@ -133,12 +145,18 @@ def alpha(val, geo): DASolver.updateDAOption() +def twist(val, geo): + for i in range(2): + geo.rot_z["bodyAxis"].coef[i] = -val[0] + + # select points pts = DVGeo.getLocalIndex(0) indexList = pts[1:4, 1, 0].flatten() PS = geo_utils.PointSelect("list", indexList) DVGeo.addLocalDV("shapey", lower=-1.0, upper=1.0, axis="y", scale=1.0, pointSelect=PS) DVGeo.addGlobalDV("alpha", [alpha0], alpha, lower=-10.0, upper=10.0, scale=1.0) +DVGeo.addGlobalDV("twist", np.zeros(1), twist, lower=-10.0, upper=10.0, scale=1.0) # DAFoam DASolver = PYDAFOAM(options=aeroOptions, comm=gcomm) @@ -219,5 +237,6 @@ def alpha(val, geo): lines = f.readlines() f.close() line4 = float(lines[4]) - if abs(line4 - 0.002349192076814971) / line4 > 1e-3: + print(line4) + if abs(line4 - 0.0022861597134965495) / line4 > 1e-3: exit(1)