diff --git a/dafoam/pyDAFoam.py b/dafoam/pyDAFoam.py index a83cfb96..6477d6ce 100755 --- a/dafoam/pyDAFoam.py +++ b/dafoam/pyDAFoam.py @@ -2275,7 +2275,7 @@ def calcTotalDerivsACT(self, objFuncName, designVarName, designVarType, dFScalin else: self.adjTotalDeriv[objFuncName][designVarName][i] = totalDerivSeq[i] - def calcTotalDerivsRegPar(self, objFuncName, designVarName, accumulateTotal=False): + def calcTotalDerivsRegPar(self, objFuncName, designVarName, dFScaling=1.0, accumulateTotal=False): xDV = self.DVGeo.getValues() nDVs = len(xDV[designVarName]) @@ -2284,34 +2284,39 @@ def calcTotalDerivsRegPar(self, objFuncName, designVarName, accumulateTotal=Fals if nDVs != nParameters: raise Error("number of parameters not valid!") - # We assume dFdRegPar is always zero - # call the total deriv - xvArray = self.vec2Array(self.xvVec) wArray = self.vec2Array(self.wVec) seedArray = self.vec2Array(self.adjVectors[objFuncName]) parameters = xDV[designVarName].copy(order="C") - productArray = np.zeros(nDVs) + totalDerivArray = np.zeros(nDVs) + dFdRegPar = np.zeros(nDVs) + + # calc dFdRegPar + self.solverAD.calcdFdRegParAD( + xvArray, wArray, parameters, objFuncName.encode(), designVarName.encode(), dFdRegPar + ) + dFdRegPar *= dFScaling # calculate dRdFieldT*Psi and save it to totalDeriv - self.solverAD.calcdRdRegParTPsiAD(xvArray, wArray, parameters, seedArray, productArray) + self.solverAD.calcdRdRegParTPsiAD(xvArray, wArray, parameters, seedArray, totalDerivArray) # all reduce because parameters is a global DV - productArray = self.comm.allreduce(productArray, op=MPI.SUM) + totalDerivArray = self.comm.allreduce(totalDerivArray, op=MPI.SUM) # totalDeriv = dFdRegPar - dRdRegParT*psi - productArray *= -1.0 + totalDerivArray = dFdRegPar - totalDerivArray # assign the total derivative to self.adjTotalDeriv if self.adjTotalDeriv[objFuncName][designVarName] is None: self.adjTotalDeriv[objFuncName][designVarName] = np.zeros(nDVs, self.dtype) - # NOTE: productArray is already in Seq + # NOTE: totalDerivArray is already in Seq because we have called all reduce in dFdRegPar + # and after calcdRdRegParTPsiAD for i in range(nDVs): if accumulateTotal is True: - self.adjTotalDeriv[objFuncName][designVarName][i] += productArray[i] + self.adjTotalDeriv[objFuncName][designVarName][i] += totalDerivArray[i] else: - self.adjTotalDeriv[objFuncName][designVarName][i] = productArray[i] + self.adjTotalDeriv[objFuncName][designVarName][i] = totalDerivArray[i] def solveAdjointUnsteady(self): """ @@ -2514,7 +2519,7 @@ def solveAdjointUnsteady(self): fieldType = designVarDict[designVarName]["fieldType"] self.calcTotalDerivsField(objFuncName, designVarName, fieldType, dFScaling, True) elif designVarDict[designVarName]["designVarType"] == "RegPar": - self.calcTotalDerivsRegPar(objFuncName, designVarName, True) + self.calcTotalDerivsRegPar(objFuncName, designVarName, dFScaling, True) else: raise Error("designVarType not valid!") @@ -4064,7 +4069,7 @@ def _initOption(self, name, value): "Expected data type is %-47s \n " "Received data type is %-47s" % (name, self.defaultOptions[name][0], type(value)) ) - + def setRegressionParameter(self, idx, val): """ Update the regression parameters diff --git a/src/adjoint/DAObjFunc/DAObjFuncVariance.C b/src/adjoint/DAObjFunc/DAObjFuncVariance.C index 423da9d4..5edd0391 100755 --- a/src/adjoint/DAObjFunc/DAObjFuncVariance.C +++ b/src/adjoint/DAObjFunc/DAObjFuncVariance.C @@ -65,6 +65,8 @@ DAObjFuncVariance::DAObjFuncVariance( objFuncDict_.readEntry("components", components_); } + timeDependentRefData_ = objFuncDict_.getLabel("timeDependentRefData"); + if (daIndex.adjStateNames.found(varName_)) { objFuncConInfo_ = {{varName_}}; @@ -75,6 +77,21 @@ DAObjFuncVariance::DAObjFuncVariance( scalar deltaT = mesh_.time().deltaT().value(); label nTimeSteps = round(endTime / deltaT); + word checkRefDataFolder; + label nRefValueInstances; + if (timeDependentRefData_) + { + // check if we can find the ref data in the endTime folder + checkRefDataFolder = Foam::name(endTime); + nRefValueInstances = nTimeSteps; + } + else + { + // check if we can find the ref data in the 0 folder + checkRefDataFolder = Foam::name(0); + nRefValueInstances = 1; + } + // check if the reference data files exist isRefData_ = 1; if (varType_ == "scalar") @@ -82,7 +99,7 @@ DAObjFuncVariance::DAObjFuncVariance( volScalarField varData( IOobject( varName_ + "Data", - Foam::name(endTime), + checkRefDataFolder, mesh_, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE), @@ -100,7 +117,7 @@ DAObjFuncVariance::DAObjFuncVariance( volVectorField varData( IOobject( varName_ + "Data", - Foam::name(endTime), + checkRefDataFolder, mesh_, IOobject::READ_IF_PRESENT, IOobject::NO_WRITE), @@ -134,15 +151,23 @@ DAObjFuncVariance::DAObjFuncVariance( { // varData file found, we need to read in the ref values for all time instances - refValue_.setSize(nTimeSteps); + refValue_.setSize(nRefValueInstances); // set refValue if (varType_ == "scalar") { - for (label n = 0; n < nTimeSteps; n++) + for (label n = 0; n < nRefValueInstances; n++) { - scalar t = (n + 1) * deltaT; - word timeName = Foam::name(t); + word timeName; + if (timeDependentRefData_) + { + scalar t = (n + 1) * deltaT; + timeName = Foam::name(t); + } + else + { + timeName = Foam::name(0); + } volScalarField varData( IOobject( @@ -210,10 +235,18 @@ DAObjFuncVariance::DAObjFuncVariance( } else if (varType_ == "vector") { - for (label n = 0; n < nTimeSteps; n++) + for (label n = 0; n < nRefValueInstances; n++) { - scalar t = (n + 1) * deltaT; - word timeName = Foam::name(t); + word timeName; + if (timeDependentRefData_) + { + scalar t = (n + 1) * deltaT; + timeName = Foam::name(t); + } + else + { + timeName = Foam::name(0); + } volVectorField varData( IOobject( @@ -339,7 +372,15 @@ void DAObjFuncVariance::calcObjFunc( const objectRegistry& db = mesh_.thisDb(); - label timeIndex = mesh_.time().timeIndex(); + label timeIndex; + if (timeDependentRefData_) + { + timeIndex = mesh_.time().timeIndex(); + } + else + { + timeIndex = 1; + } if (varName_ == "wallShearStress") { diff --git a/src/adjoint/DAObjFunc/DAObjFuncVariance.H b/src/adjoint/DAObjFunc/DAObjFuncVariance.H index a09ebb83..a9dbff9c 100755 --- a/src/adjoint/DAObjFunc/DAObjFuncVariance.H +++ b/src/adjoint/DAObjFunc/DAObjFuncVariance.H @@ -61,6 +61,9 @@ protected: /// whether we find the reference data label isRefData_; + /// whether the ref data is time dependent if yes we need data in all time folders otherwise get it from the 0 folder + label timeDependentRefData_; + /// DATurbulenceModel object const DATurbulenceModel& daTurb_; diff --git a/src/adjoint/DASolver/DASolver.C b/src/adjoint/DASolver/DASolver.C index 5853a2cc..685dabe8 100644 --- a/src/adjoint/DASolver/DASolver.C +++ b/src/adjoint/DASolver/DASolver.C @@ -5305,6 +5305,10 @@ void DASolver::updateOFField(const Vec wVec) } this->setPrimalBoundaryConditions(printInfo); daFieldPtr_->stateVec2OFField(wVec); + + // if we have regression models, we also need to update them because they will update the fields + this->regressionModelCompute(); + // We need to call correctBC multiple times to reproduce // the exact residual, this is needed for some boundary conditions // and intermediate variables (e.g., U for inletOutlet, nut with wall functions) @@ -5328,6 +5332,10 @@ void DASolver::updateOFField(const scalar* states) } this->setPrimalBoundaryConditions(printInfo); daFieldPtr_->state2OFField(states); + + // if we have regression models, we also need to update them because they will update the fields + this->regressionModelCompute(); + // We need to call correctBC multiple times to reproduce // the exact residual, this is needed for some boundary conditions // and intermediate variables (e.g., U for inletOutlet, nut with wall functions) @@ -5758,6 +5766,158 @@ void DASolver::calcdFdXvAD( #endif } +void DASolver::calcdFdRegParAD( + const double* volCoords, + const double* states, + const double* parameters, + const word objFuncName, + const word designVarName, + double* dFdRegPar) +{ +#ifdef CODI_AD_REVERSE + /* + Description: + Compute dFdRegPar using reverse-mode AD + + Input: + + xvVec: the volume mesh coordinate vector + + wVec: the state variable vector + + objFuncName: the name of the objective function + + designVarName: name of the design variable + + Output: + dFdRegPar: dF/dRegPar + */ + + Info << "Calculating dFdRegPar using reverse-mode AD" << endl; + + scalar* volCoordsArray = new scalar[daIndexPtr_->nLocalXv]; + for (label i = 0; i < daIndexPtr_->nLocalXv; i++) + { + volCoordsArray[i] = volCoords[i]; + } + + scalar* statesArray = new scalar[daIndexPtr_->nLocalAdjointStates]; + for (label i = 0; i < daIndexPtr_->nLocalAdjointStates; i++) + { + statesArray[i] = states[i]; + } + + label nParameters = this->getNRegressionParameters(); + scalar* parametersArray = new scalar[nParameters]; + for (label i = 0; i < nParameters; i++) + { + parametersArray[i] = parameters[i]; + // NOTE: we also zero out the dFdRegPar array + dFdRegPar[i] = 0; + } + // for each objFunc part + scalarList dFdRegParPart(nParameters, 0.0); + + this->updateOFMesh(volCoordsArray); + this->updateOFField(statesArray); + + // get the subDict for this objective function + dictionary objFuncSubDict = + daOptionPtr_->getAllOptions().subDict("objFunc").subDict(objFuncName); + + // loop over all parts for this objFuncName + forAll(objFuncSubDict.toc(), idxJ) + { + // get the subDict for this part + word objFuncPart = objFuncSubDict.toc()[idxJ]; + + // get objFunc from daObjFuncPtrList_ + label objIndx = this->getObjFuncListIndex(objFuncName, objFuncPart); + DAObjFunc& daObjFunc = daObjFuncPtrList_[objIndx]; + + // reset tape + this->globalADTape_.reset(); + // activate tape, start recording + this->globalADTape_.setActive(); + // register inputs + for (label i = 0; i < nParameters; i++) + { + this->globalADTape_.registerInput(parametersArray[i]); + } + for (label i = 0; i < nParameters; i++) + { + this->setRegressionParameter(i, parametersArray[i]); + } + + // update the BC, this func will also call regressionModelCompute + this->updateStateBoundaryConditions(); + + // compute the objective function + scalar fRef = daObjFunc.getObjFuncValue(); + // register f as the output + this->globalADTape_.registerOutput(fRef); + // stop recording + this->globalADTape_.setPassive(); + + // Note: since we used reduced objFunc, we only need to + // assign the seed for master proc + if (Pstream::master()) + { + fRef.setGradient(1.0); + } + // evaluate tape to compute derivative + this->globalADTape_.evaluate(); + + // get the partials from the inputs. + for (label i = 0; i < nParameters; i++) + { + dFdRegParPart[i] = parametersArray[i].getGradient(); + } + + for (label i = 0; i < nParameters; i++) + { + // we need to reduce all the contributions from processors because regPar is a serial input + reduce(dFdRegParPart[i], sumOp()); + + // we need to add dFd*Part to dFd* because we want to sum + // all dFd*Part for all parts of this objFuncName. + dFdRegPar[i] += dFdRegParPart[i].getValue(); + } + + // need to clear adjoint and tape after the computation is done! + this->globalADTape_.clearAdjoints(); + this->globalADTape_.reset(); + + // ********************************************************************************************** + // clean up OF vars's AD seeds by deactivating the inputs and call the forward func one more time + // ********************************************************************************************** + + for (label i = 0; i < nParameters; i++) + { + this->globalADTape_.deactivateValue(parametersArray[i]); + } + for (label i = 0; i < nParameters; i++) + { + this->setRegressionParameter(i, parametersArray[i]); + } + + this->updateStateBoundaryConditions(); + fRef = daObjFunc.getObjFuncValue(); + + if (daOptionPtr_->getOption