Skip to content

Commit

Permalink
Fixed a bug to allow mulitple parts for unsteady objs (#536)
Browse files Browse the repository at this point in the history 
* Fixed a bug to allow mulitple parts for unsteady objs.

* Updated test scripts.
  • Loading branch information
@friedenhe
friedenhe authored Dec 7, 2023
1 parent 227bad5 commit 3559006
Show file tree
Hide file tree
Showing 27 changed files with 102 additions and 105 deletions.
3 changes: 2 additions & 1 deletion dafoam/pyDAFoam.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -465,6 +465,7 @@ def __init__(self):
"periodicity": -1.0,
"objFuncStartTime": -1.0,
"objFuncEndTime": -1.0,
"objFuncTimeOperator": "None",
"PCMatPrecomputeInterval": 100,
"PCMatUpdateInterval": 1,
"reduceIO": True,
Expand Down Expand Up @@ -2415,7 +2416,7 @@ def solveAdjointUnsteady(self):
# on the state at this time index
dFScaling = 1.0
if n >= objFuncStartTimeIndex and n <= objFuncEndTimeIndex:
dFScaling = self.solver.getObjFuncUnsteadyScaling(objFuncName.encode())
dFScaling = self.solver.getObjFuncUnsteadyScaling()
else:
dFScaling = 0.0

Expand Down
8 changes: 0 additions & 8 deletions src/adjoint/DAObjFunc/DAObjFunc.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,9 +92,6 @@ protected:
/// the connectivity information for the objective function
List<List<word>> objFuncConInfo_;

/// the unsteady objective function type like sum or average
word timeOperator_ = "None";

public:
/// Runtime type information
TypeName("DAObjFunc");
Expand Down Expand Up @@ -234,11 +231,6 @@ public:
return objFuncConInfo_;
}

const word& getObjFuncTimeOperator() const
{
return timeOperator_;
}

/// expSumKS stores sum[exp(coeffKS*x_i)] for KS function which will be used to scale dFdW
scalar expSumKS = 1.0;

Expand Down
2 changes: 0 additions & 2 deletions src/adjoint/DAObjFunc/DAObjFuncCenterOfPressure.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -71,8 +71,6 @@ DAObjFuncCenterOfPressure::DAObjFuncCenterOfPressure(
center_[2] = centerList[2];

objFuncDict_.readEntry<scalar>("scale", scale_);

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
2 changes: 0 additions & 2 deletions src/adjoint/DAObjFunc/DAObjFuncFieldInversion.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,8 +74,6 @@ DAObjFuncFieldInversion::DAObjFuncFieldInversion(
{
objFuncConInfo_ = {{}}; // level 0
}

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
1 change: 0 additions & 1 deletion src/adjoint/DAObjFunc/DAObjFuncForce.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -115,7 +115,6 @@ DAObjFuncForce::DAObjFuncForce(
// now replace nut with the corrected name for the selected turbulence model
daModel.correctModelStates(objFuncConInfo_[0]);

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
1 change: 0 additions & 1 deletion src/adjoint/DAObjFunc/DAObjFuncLocation.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -110,7 +110,6 @@ DAObjFuncLocation::DAObjFuncLocation(
}
}

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
1 change: 0 additions & 1 deletion src/adjoint/DAObjFunc/DAObjFuncMass.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,6 @@ DAObjFuncMass::DAObjFuncMass(

objFuncDict_.readEntry<scalar>("scale", scale_);

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
1 change: 0 additions & 1 deletion src/adjoint/DAObjFunc/DAObjFuncMassFlowRate.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -57,7 +57,6 @@ DAObjFuncMassFlowRate::DAObjFuncMassFlowRate(

objFuncDict_.readEntry<scalar>("scale", scale_);

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
1 change: 0 additions & 1 deletion src/adjoint/DAObjFunc/DAObjFuncMeshQualityKS.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,6 @@ DAObjFuncMeshQualityKS::DAObjFuncMeshQualityKS(
Info << "includeFaceList " << includeFaceList_ << endl;
}

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
1 change: 0 additions & 1 deletion src/adjoint/DAObjFunc/DAObjFuncMoment.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,6 @@ DAObjFuncMoment::DAObjFuncMoment(
// now replace nut with the corrected name for the selected turbulence model
daModel.correctModelStates(objFuncConInfo_[0]);

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
2 changes: 0 additions & 2 deletions src/adjoint/DAObjFunc/DAObjFuncPatchMean.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -45,8 +45,6 @@ DAObjFuncPatchMean::DAObjFuncPatchMean(
objFuncDict_.readEntry<word>("varType", varType_);

objFuncDict_.readEntry<label>("component", component_);

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
2 changes: 0 additions & 2 deletions src/adjoint/DAObjFunc/DAObjFuncPower.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -90,8 +90,6 @@ DAObjFuncPower::DAObjFuncPower(

// now replace nut with the corrected name for the selected turbulence model
daModel.correctModelStates(objFuncConInfo_[0]);

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
2 changes: 0 additions & 2 deletions src/adjoint/DAObjFunc/DAObjFuncTotalPressure.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -56,8 +56,6 @@ DAObjFuncTotalPressure::DAObjFuncTotalPressure(
#endif

objFuncDict_.readEntry<scalar>("scale", scale_);

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
1 change: 0 additions & 1 deletion src/adjoint/DAObjFunc/DAObjFuncTotalPressureRatio.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,7 +65,6 @@ DAObjFuncTotalPressureRatio::DAObjFuncTotalPressureRatio(
Info << "gamma " << gamma_ << endl;
}

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
1 change: 0 additions & 1 deletion src/adjoint/DAObjFunc/DAObjFuncTotalTemperatureRatio.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,6 @@ DAObjFuncTotalTemperatureRatio::DAObjFuncTotalTemperatureRatio(
Info << "gamma " << gamma_ << endl;
}

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
2 changes: 0 additions & 2 deletions src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -51,8 +51,6 @@ DAObjFuncVariableVolSum::DAObjFuncVariableVolSum(

objFuncDict_.readEntry<label>("divByTotalVol", divByTotalVol_);

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");

if (daIndex.adjStateNames.found(varName_))
{
objFuncConInfo_ = {{varName_}};
Expand Down
2 changes: 0 additions & 2 deletions src/adjoint/DAObjFunc/DAObjFuncVariance.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -65,8 +65,6 @@ DAObjFuncVariance::DAObjFuncVariance(
objFuncDict_.readEntry<labelList>("components", components_);
}

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");

if (daIndex.adjStateNames.found(varName_))
{
objFuncConInfo_ = {{varName_}};
Expand Down
1 change: 0 additions & 1 deletion src/adjoint/DAObjFunc/DAObjFuncVonMisesStressKS.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -47,7 +47,6 @@ DAObjFuncVonMisesStressKS::DAObjFuncVonMisesStressKS(

objFuncDict_.readEntry<scalar>("coeffKS", coeffKS_);

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");
}

/// calculate the value of objective function
Expand Down
2 changes: 0 additions & 2 deletions src/adjoint/DAObjFunc/DAObjFuncWallHeatFlux.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,8 +79,6 @@ DAObjFuncWallHeatFlux::DAObjFuncWallHeatFlux(

objFuncDict_.readEntry<scalar>("scale", scale_);

timeOperator_ = objFuncDict.lookupOrDefault<word>("timeOperator", "None");

#ifdef CompressibleFlow

// setup the connectivity for heat flux, this is needed in Foam::DAJacCondFdW
Expand Down
84 changes: 66 additions & 18 deletions src/adjoint/DASolver/DASolver.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -173,26 +173,30 @@ void DASolver::calcUnsteadyObjFuncs()

scalar objFuncValue = 0.0;

word timeOperator = daOptionPtr_->getSubDictOption<word>("unsteadyAdjoint", "objFuncTimeOperator");

forAll(daObjFuncPtrList_, idxI)
{
DAObjFunc& daObjFunc = daObjFuncPtrList_[idxI];
word objFuncName = daObjFunc.getObjFuncName();
if (daObjFunc.getObjFuncTimeOperator() == "None")
word objFuncPart = daObjFunc.getObjFuncPart();
word uKey = objFuncName + objFuncPart;
if (timeOperator == "None")
{
FatalErrorIn("") << "calcUnsteadyObjFuncs is called but the timeOperator is not set!!! Options are: average or sum"
<< abort(FatalError);
}
else if (daObjFunc.getObjFuncTimeOperator() == "sum")
else if (timeOperator == "sum")
{
label startTimeIndex = this->getUnsteadyObjFuncStartTimeIndex();
label endTimeIndex = this->getUnsteadyObjFuncEndTimeIndex();
label timeIndex = runTimePtr_->timeIndex();
if (timeIndex >= startTimeIndex && timeIndex <= endTimeIndex)
{
unsteadyObjFuncs_[objFuncName] += this->getObjFuncValue(objFuncName);
unsteadyObjFuncs_[uKey] += daObjFunc.getObjFuncValue();
}
}
else if (daObjFunc.getObjFuncTimeOperator() == "average")
else if (timeOperator == "average")
{
// calculate the average on the fly, i.e., moving average
label startTimeIndex = this->getUnsteadyObjFuncStartTimeIndex();
Expand All @@ -201,14 +205,14 @@ void DASolver::calcUnsteadyObjFuncs()
if (timeIndex >= startTimeIndex && timeIndex <= endTimeIndex)
{
label n = timeIndex - startTimeIndex + 1;
scalar objFuncVal = this->getObjFuncValue(objFuncName);
unsteadyObjFuncs_[objFuncName] = (unsteadyObjFuncs_[objFuncName] * (n - 1) + objFuncVal) / n;
scalar objFuncVal = daObjFunc.getObjFuncValue();
unsteadyObjFuncs_[uKey] = (unsteadyObjFuncs_[uKey] * (n - 1) + objFuncVal) / n;
}
}
else
{
FatalErrorIn("calcUnsteadyObjFuncs") << "timeOperator not valid! Options are None, average, or sum"
<< abort(FatalError);
FatalErrorIn("") << "calcUnsteadyObjFuncs is called but the timeOperator is not set!!! Options are: average or sum"
<< abort(FatalError);
}
}
}
Expand All @@ -228,18 +232,22 @@ void DASolver::printAllObjFuncs()
<< abort(FatalError);
}

word timeOperator = daOptionPtr_->getSubDictOption<word>("unsteadyAdjoint", "objFuncTimeOperator");

forAll(daObjFuncPtrList_, idxI)
{
DAObjFunc& daObjFunc = daObjFuncPtrList_[idxI];
word objFuncName = daObjFunc.getObjFuncName();
word objFuncPart = daObjFunc.getObjFuncPart();
word uKey = objFuncName + objFuncPart;
scalar objFuncVal = daObjFunc.getObjFuncValue();
Info << objFuncName
<< "-" << daObjFunc.getObjFuncPart()
<< "-" << objFuncPart
<< "-" << daObjFunc.getObjFuncType()
<< ": " << objFuncVal;
if (daObjFunc.getObjFuncTimeOperator() == "average" || daObjFunc.getObjFuncTimeOperator() == "sum")
if (timeOperator == "average" || timeOperator == "sum")
{
Info << " Unsteady " << daObjFunc.getObjFuncTimeOperator() << " " << unsteadyObjFuncs_[objFuncName];
Info << " Unsteady " << timeOperator << " " << unsteadyObjFuncs_[uKey];
}
#ifdef CODI_AD_FORWARD

Expand All @@ -254,14 +262,52 @@ void DASolver::printAllObjFuncs()

if (daOptionPtr_->getSubDictOption<word>("unsteadyAdjoint", "mode") == "timeAccurate")
{
Info << " Unsteady Deriv: " << unsteadyObjFuncs_[objFuncName].getGradient();
Info << " Unsteady Deriv: " << unsteadyObjFuncs_[uKey].getGradient();
}
}
#endif
Info << endl;
}
}

scalar DASolver::getObjFuncValueUnsteady(const word objFuncName)
{
/*
Description:
Return the value of the objective function for unsteady cases
NOTE: we will sum up all the parts in objFuncName
Input:
objFuncName: the name of the objective function
Output:
objFuncValue: the value of the objective
*/
if (daObjFuncPtrList_.size() == 0)
{
FatalErrorIn("printAllObjFuncs") << "daObjFuncPtrList_.size() ==0... "
<< "Forgot to call setDAObjFuncList?"
<< abort(FatalError);
}

scalar objFuncValue = 0.0;

forAll(daObjFuncPtrList_, idxI)
{
DAObjFunc& daObjFunc = daObjFuncPtrList_[idxI];
word objFuncNameI = daObjFunc.getObjFuncName();
word objFuncPart = daObjFunc.getObjFuncPart();
word uKey = objFuncNameI + objFuncPart;

if (objFuncNameI == objFuncName)
{
objFuncValue += unsteadyObjFuncs_[uKey];
}
}

return objFuncValue;
}

scalar DASolver::getObjFuncValue(const word objFuncName)
{
/*
Expand Down Expand Up @@ -398,24 +444,26 @@ void DASolver::setDAObjFuncList()
{
DAObjFunc& daObjFunc = daObjFuncPtrList_[idxI];
word objFuncName = daObjFunc.getObjFuncName();
unsteadyObjFuncs_.set(objFuncName, 0.0);
word objFuncPart = daObjFunc.getObjFuncPart();
word uKey = objFuncName + objFuncPart;
unsteadyObjFuncs_.set(uKey, 0.0);

word timeOperator = daObjFunc.getObjFuncTimeOperator();
word timeOperator = daOptionPtr_->getSubDictOption<word>("unsteadyAdjoint", "objFuncTimeOperator");
if (timeOperator == "None" || timeOperator == "sum")
{
unsteadyObjFuncsScaling_.set(objFuncName, 1.0);
unsteadyObjFuncsScaling_ = 1.0;
}
else if (timeOperator == "average")
{
label startTimeIndex = this->getUnsteadyObjFuncStartTimeIndex();
label endTimeIndex = this->getUnsteadyObjFuncEndTimeIndex();
label nInstances = endTimeIndex - startTimeIndex + 1;
unsteadyObjFuncsScaling_.set(objFuncName, 1.0 / nInstances);
unsteadyObjFuncsScaling_ = 1.0 / nInstances;
}
else
{
FatalErrorIn("calcUnsteadyObjFuncs") << "timeOperator not valid! Options are None, average, or sum"
<< abort(FatalError);
FatalErrorIn("setDAObjFuncList") << "timeOperator not valid! Options are None, average, or sum"
<< abort(FatalError);
}
}
}
Expand Down
15 changes: 7 additions & 8 deletions src/adjoint/DASolver/DASolver.H
View file
Original file line number Diff line number Diff line change
Expand Up @@ -207,7 +207,7 @@ protected:
HashTable<scalar> unsteadyObjFuncs_;

/// a scaling factor related to the unsteady obj func timeOperator e.g., in dFdW calculation
HashTable<scalar> unsteadyObjFuncsScaling_;
scalar unsteadyObjFuncsScaling_;

public:
/// Runtime type information
Expand Down Expand Up @@ -780,23 +780,22 @@ public:
{
DAObjFunc& daObjFunc = daObjFuncPtrList_[idxI];
word objFuncName = daObjFunc.getObjFuncName();
unsteadyObjFuncs_[objFuncName] = 0.0;
word objFuncPart = daObjFunc.getObjFuncPart();
word uKey = objFuncName + objFuncPart;
unsteadyObjFuncs_[uKey] = 0.0;
}
}

/// calculate and save the unsteady objective function
void calcUnsteadyObjFuncs();

/// return the unsteady objFunc values
scalar getObjFuncValueUnsteady(const word objFuncName)
{
return unsteadyObjFuncs_[objFuncName];
}
scalar getObjFuncValueUnsteady(const word objFuncName);

/// get the objective func scaling for unsteady adjoint
scalar getObjFuncUnsteadyScaling(const word objFuncName)
scalar getObjFuncUnsteadyScaling()
{
return unsteadyObjFuncsScaling_[objFuncName];
return unsteadyObjFuncsScaling_;
}

/// get the unsteady adjoint start time index
Expand Down
Loading

0 comments on commit 3559006

Please sign in to comment.