Skip to content

Commit

Permalink
Added an option for maxRadius location objFunc (#509)
Browse files Browse the repository at this point in the history 
* Added an option for maxRadius.

* Updated script.
  • Loading branch information
@friedenhe
friedenhe authored Oct 28, 2023
1 parent 3956592 commit dfc14a0
Show file tree
Hide file tree
Showing 4 changed files with 135 additions and 24 deletions.
110 changes: 88 additions & 22 deletions src/adjoint/DAObjFunc/DAObjFuncLocation.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<scalar>());

if (fabs(maxRGlobal - maxR) < 1e-8)
{
maxRPatchI_ = maxRPatchI;
maxRFaceI_ = maxRFaceI;
}
}
}

/// calculate the value of objective function
Expand Down Expand Up @@ -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();
Expand All @@ -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<scalar>());

// 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<scalar>());
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<scalar>());

objFuncValue = log(objValTmp) / coeffKS_;
objFuncValue = radius;
}
else
{
FatalErrorIn("DAObjFuncLocation") << "mode: " << mode_ << " not supported!"
<< "Options are: maxRadius"
<< abort(FatalError);
}

return;
}
Expand Down
7 changes: 7 additions & 0 deletions src/adjoint/DAObjFunc/DAObjFuncLocation.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
19 changes: 19 additions & 0 deletions tests/refs/DAFoam_Test_DASimpleFoamADRef.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -17,27 +19,44 @@ 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
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
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
Expand Down
23 changes: 21 additions & 2 deletions tests/runTests_DASimpleFoamAD.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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},
Expand All @@ -107,6 +118,7 @@
"designVar": {
"shapey": {"designVarType": "FFD"},
"alpha": {"designVarType": "AOA", "patches": ["inout"], "flowAxis": "x", "normalAxis": "y"},
"twist": {"designVarType": "FFD"}
},
}

Expand All @@ -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)
Expand Down Expand Up @@ -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)

0 comments on commit dfc14a0

Please sign in to comment.