From 172868398d64abfde56aa82c2a9e72f387561dd4 Mon Sep 17 00:00:00 2001 From: Ping He Date: Wed, 27 Dec 2023 21:26:34 -0600 Subject: [PATCH 1/8] Added DARegression and its derivs. --- dafoam/pyDAFoam.py | 81 ++++++- src/adjoint/DAModel/DAModel.C | 22 ++ src/adjoint/DAModel/DAModel.H | 6 +- .../DATurbulenceModel/DASpalartAllmaras.C | 21 ++ .../DATurbulenceModel/DASpalartAllmaras.H | 3 + .../DATurbulenceModel/DASpalartAllmarasFv3.C | 25 +- .../DATurbulenceModel/DASpalartAllmarasFv3.H | 3 + .../DATurbulenceModel/DATurbulenceModel.C | 14 +- .../DATurbulenceModel/DATurbulenceModel.H | 5 +- src/adjoint/DARegression/DARegression.C | 220 ++++++++++++++++++ src/adjoint/DARegression/DARegression.H | 110 +++++++++ .../DASolver/DAPimpleFoam/DAPimpleFoam.C | 18 ++ .../DASolver/DAPimpleFoam/DAPimpleFoam.H | 3 + src/adjoint/DASolver/DASolver.C | 122 +++++++++- src/adjoint/DASolver/DASolver.H | 36 +++ src/adjoint/Make/files_Incompressible | 2 + src/include/createAdjointIncompressible.H | 7 + src/include/setForwardADSeeds.H | 13 ++ src/pyDASolvers/DASolvers.H | 30 +++ src/pyDASolvers/pyDASolvers.pyx | 39 ++++ 20 files changed, 773 insertions(+), 7 deletions(-) create mode 100755 src/adjoint/DARegression/DARegression.C create mode 100755 src/adjoint/DARegression/DARegression.H diff --git a/dafoam/pyDAFoam.py b/dafoam/pyDAFoam.py index 0f1c2873..a83cfb96 100755 --- a/dafoam/pyDAFoam.py +++ b/dafoam/pyDAFoam.py @@ -508,6 +508,17 @@ def __init__(self): ## constrainHbyA back to the primal and adjoint solvers. self.useConstrainHbyA = False + ## parameters for regression models + self.regressionModel = { + "active": False, + "modelType": "neuralNetwork", + "inputNames": ["None"], + "outputName": "None", + "hiddenLayerNeurons": [0], + "outputShift": 0.0, + "outputScale": 1.0, + } + # ********************************************************************************************* # ************************************ Advance Options **************************************** # ********************************************************************************************* @@ -2264,6 +2275,44 @@ def calcTotalDerivsACT(self, objFuncName, designVarName, designVarType, dFScalin else: self.adjTotalDeriv[objFuncName][designVarName][i] = totalDerivSeq[i] + def calcTotalDerivsRegPar(self, objFuncName, designVarName, accumulateTotal=False): + + xDV = self.DVGeo.getValues() + nDVs = len(xDV[designVarName]) + + nParameters = self.solver.getNRegressionParameters() + 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) + + # calculate dRdFieldT*Psi and save it to totalDeriv + self.solverAD.calcdRdRegParTPsiAD(xvArray, wArray, parameters, seedArray, productArray) + # all reduce because parameters is a global DV + productArray = self.comm.allreduce(productArray, op=MPI.SUM) + + # totalDeriv = dFdRegPar - dRdRegParT*psi + productArray *= -1.0 + + # 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 + + for i in range(nDVs): + if accumulateTotal is True: + self.adjTotalDeriv[objFuncName][designVarName][i] += productArray[i] + else: + self.adjTotalDeriv[objFuncName][designVarName][i] = productArray[i] + def solveAdjointUnsteady(self): """ Run adjoint solver to compute the total derivs for unsteady solvers @@ -2274,6 +2323,8 @@ def solveAdjointUnsteady(self): wVec: vector that contains all the state variables + designVariable: a dictionary that contain the design variable values + Output: ------- self.adjTotalDeriv: the dict contains all the total derivative vectors @@ -2462,6 +2513,8 @@ def solveAdjointUnsteady(self): elif designVarDict[designVarName]["designVarType"] == "Field": fieldType = designVarDict[designVarName]["fieldType"] self.calcTotalDerivsField(objFuncName, designVarName, fieldType, dFScaling, True) + elif designVarDict[designVarName]["designVarType"] == "RegPar": + self.calcTotalDerivsRegPar(objFuncName, designVarName, True) else: raise Error("designVarType not valid!") @@ -2751,6 +2804,15 @@ def solveAdjoint(self): self.calcTotalDerivsField(objFuncName, designVarName, fieldType) else: raise Error("For Field design variable type, we only support useAD->mode=reverse") + ################### RegPar: regression model's parameters ################### + elif designVarDict[designVarName]["designVarType"] == "RegPar": + if self.getOption("useAD")["mode"] == "reverse": + # loop over all objectives + for objFuncName in objFuncDict: + if objFuncName in self.objFuncNames4Adj: + self.calcTotalDerivsRegPar(objFuncName, designVarName) + else: + raise Error("For Field design variable type, we only support useAD->mode=reverse") else: raise Error("designVarType %s not supported!" % designVarDict[designVarName]["designVarType"]) @@ -3151,7 +3213,7 @@ def _initSolver(self): adjMode = self.getOption("unsteadyAdjoint")["mode"] if adjMode == "hybrid": self.initTimeInstanceMats() - + Info("Init solver done! ElapsedClockTime %f s" % self.solver.getElapsedClockTime()) Info("Init solver done!. ElapsedCpuTime %f s" % self.solver.getElapsedCpuTime()) @@ -4002,6 +4064,23 @@ 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 + + Parameters + ---------- + idx: int + the index of the parameter + + val, float + the parameter value to set + """ + + self.solver.setRegressionParameter(idx, val) + if self.getOption("useAD")["mode"] in ["forward", "reverse"]: + self.solverAD.setRegressionParameter(idx, val) def setFieldValue4GlobalCellI(self, fieldName, val, globalCellI, compI=0): """ diff --git a/src/adjoint/DAModel/DAModel.C b/src/adjoint/DAModel/DAModel.C index 4fe38cf1..91c93cd9 100755 --- a/src/adjoint/DAModel/DAModel.C +++ b/src/adjoint/DAModel/DAModel.C @@ -292,6 +292,28 @@ void DAModel::getTurbProdTerm(scalarList& prodTerm) const } +void DAModel::getTurbProdOverDestruct(scalarList& PoD) const +{ + /* + Description: + Return the value of the production/destruction term from the turbulence model + */ + +#ifndef SolidDASolver + if (hasTurbulenceModel_) + { + DATurbulenceModel& daTurb = const_cast( + mesh_.thisDb().lookupObject("DATurbulenceModel")); + daTurb.getTurbProdOverDestruct(PoD); + } + + if (hasRadiationModel_) + { + } +#endif + +} + #ifdef CompressibleFlow const fluidThermo& DAModel::getThermo() const { diff --git a/src/adjoint/DAModel/DAModel.H b/src/adjoint/DAModel/DAModel.H index 5c074378..ba2b610e 100755 --- a/src/adjoint/DAModel/DAModel.H +++ b/src/adjoint/DAModel/DAModel.H @@ -105,9 +105,12 @@ public: /// update intermediate variables that are dependent on the model states void updateIntermediateVariables(); - /// return the value of the production term from the turbulence model + /// return the value of the production term from the turbulence model void getTurbProdTerm(scalarList& prodTerm) const; + /// return the value of the destruction term from the turbulence model + void getTurbProdOverDestruct(scalarList& PoD) const; + #ifndef SolidDASolver /// get a reference to DATurbulenceModel const DATurbulenceModel& getDATurbulenceModel() const; @@ -117,7 +120,6 @@ public: /// return thermo_ from DATurbulenceModel const fluidThermo& getThermo() const; #endif - }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.C b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.C index 977ba162..16de8cec 100755 --- a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.C +++ b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.C @@ -518,6 +518,27 @@ void DASpalartAllmaras::getFvMatrixFields( upper = nuTildaEqn.upper(); lower = nuTildaEqn.lower(); } + +void DASpalartAllmaras::getTurbProdOverDestruct(scalarList& PoD) const +{ + /* + Description: + Return the value of the production over destruction term from the turbulence model + */ + + const volScalarField chi(this->chi()); + const volScalarField fv1(this->fv1(chi)); + + const volScalarField Stilda(this->Stilda(chi, fv1)); + + volScalarField P = Cb1_ * phase_ * rho_ * Stilda * nuTilda_; + volScalarField D = Cw1_ * phase_ * rho_ * fw(Stilda) * sqr(nuTilda_ / y_); + + forAll(P, cellI) + { + PoD[cellI] = P[cellI] / D[cellI]; + } +} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.H b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.H index befc48f1..1f1b485a 100755 --- a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.H +++ b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmaras.H @@ -146,6 +146,9 @@ public: scalarField& diag, scalarField& upper, scalarField& lower); + + /// return the value of the destruction term from the turbulence model + virtual void getTurbProdOverDestruct(scalarList& PoD) const; }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.C b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.C index 97a6edc9..8382c61e 100755 --- a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.C +++ b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.C @@ -737,13 +737,36 @@ void DASpalartAllmarasFv3::getFvMatrixFields( - Cb2_ / sigmaNut_ * phase_ * rho_ * magSqr(fvc::grad(nuTilda_)) == Cb1_ * phase_ * rho_ * Stilda * nuTilda_ * betaFI_ - fvm::Sp(Cw1_ * phase_ * rho_ * fw(Stilda) * nuTilda_ / sqr(y_), nuTilda_)); - + nuTildaEqn.relax(); diag = nuTildaEqn.D(); upper = nuTildaEqn.upper(); lower = nuTildaEqn.lower(); } + +void DASpalartAllmarasFv3::getTurbProdOverDestruct(scalarList& PoD) const +{ + /* + Description: + Return the value of the production over destruction term from the turbulence model + */ + + const volScalarField chi(this->chi()); + const volScalarField fv1(this->fv1(chi)); + + const volScalarField Stilda( + this->fv3(chi, fv1) * ::sqrt(2.0) * mag(skew(fvc::grad(U_))) + + this->fv2(chi, fv1) * nuTilda_ / sqr(kappa_ * y_)); + + volScalarField P = Cb1_ * phase_ * rho_ * Stilda * nuTilda_; + volScalarField D = Cw1_ * phase_ * rho_ * fw(Stilda) * sqr(nuTilda_ / y_); + + forAll(P, cellI) + { + PoD[cellI] = P[cellI] / D[cellI]; + } +} // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam diff --git a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.H b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.H index 91d2fef4..09702d98 100755 --- a/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.H +++ b/src/adjoint/DAModel/DATurbulenceModel/DASpalartAllmarasFv3.H @@ -168,6 +168,9 @@ public: scalarField& diag, scalarField& upper, scalarField& lower); + + /// return the value of the destruction term from the turbulence model + virtual void getTurbProdOverDestruct(scalarList& PoD) const; }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.C b/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.C index 2ed04794..71f81b48 100755 --- a/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.C +++ b/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.C @@ -487,7 +487,19 @@ void DATurbulenceModel::getTurbProdTerm(scalarList& prodTerm) const Return the value of the production term from the turbulence model */ - FatalErrorIn("DATurbulenceModel::getSAProdTerm") + FatalErrorIn("DATurbulenceModel::getTurbProdTerm") + << "Child class not implemented!" + << abort(FatalError); +} + +void DATurbulenceModel::getTurbProdOverDestruct(scalarList& PoD) const +{ + /* + Description: + Return the value of the production over destruction term from the turbulence model + */ + + FatalErrorIn("DATurbulenceModel::getTurbProdOverDestruct") << "Child class not implemented!" << abort(FatalError); } diff --git a/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.H b/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.H index dc936934..e08bb15c 100755 --- a/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.H +++ b/src/adjoint/DAModel/DATurbulenceModel/DATurbulenceModel.H @@ -211,9 +211,12 @@ public: /// solve the residual equations and update the state virtual void correct(label printToScreen) = 0; - /// return the value of the production term from the Spalart Allmaras model + /// return the value of the production term from the turbulence model virtual void getTurbProdTerm(scalarList& prodTerm) const; + /// return the value of the destruction term from the turbulence model + virtual void getTurbProdOverDestruct(scalarList& PoD) const; + /// dev terms tmp devRhoReff() const; diff --git a/src/adjoint/DARegression/DARegression.C b/src/adjoint/DARegression/DARegression.C new file mode 100755 index 00000000..f77c8391 --- /dev/null +++ b/src/adjoint/DARegression/DARegression.C @@ -0,0 +1,220 @@ +/*---------------------------------------------------------------------------*\ + + DAFoam : Discrete Adjoint with OpenFOAM + Version : v3 + +\*---------------------------------------------------------------------------*/ + +#include "DARegression.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +DARegression::DARegression( + const fvMesh& mesh, + const DAOption& daOption, + const DAModel& daModel) + : mesh_(mesh), + daOption_(daOption), + daModel_(daModel) +{ + dictionary regSubDict = daOption.getAllOptions().subDict("regressionModel"); + regSubDict.readEntry("modelType", modelType_); + regSubDict.readEntry("inputNames", inputNames_); + regSubDict.readEntry("outputName", outputName_); + + regSubDict.readEntry("outputShift", outputShift_); + + regSubDict.readEntry("outputScale", outputScale_); + + if (modelType_ == "neuralNetwork") + { + regSubDict.readEntry("hiddenLayerNeurons", hiddenLayerNeurons_); + } + + // initialize parameters and give it large values + label nParameters = this->nParameters(); + parameters_.setSize(nParameters); + forAll(parameters_, idxI) + { + parameters_[idxI] = 1e16; + } +} + +void DARegression::compute() +{ + /* + Description: + Calculate the prescribed output field based on the prescribed input fields using a regression model. + We support only neural network at this moment + + Input: + parameters_: + the parameters for the regression. For the neural network model, these parameters + are the weights and biases. We order them in an 1D array, starting from the first input's weight and biases. + + daOption Dict: regressionModel + inputNames: a list of volScalarFields prescribed by inputFields + hiddenLayerNeurons: number of neurons for each hidden layer. + example: {5, 3, 4} means three hidden layers with 5, 3, and 4 neurons. + + Output: + a volScalarField prescribed by outputName + */ + + volScalarField& outputField = const_cast(mesh_.thisDb().lookupObject(outputName_)); + + if (modelType_ == "neuralNetwork") + { + label nHiddenLayers = hiddenLayerNeurons_.size(); + + List> inputFields; + inputFields.setSize(inputNames_.size()); + + forAll(inputNames_, idxI) + { + inputFields[idxI].setSize(mesh_.nCells()); + word inputName = inputNames_[idxI]; + if (inputName == "SoQ") + { + const volVectorField& U = mesh_.thisDb().lookupObject("U"); + const tmp tgradU(fvc::grad(U)); + const volTensorField& gradU = tgradU(); + volScalarField magOmega = mag(skew(gradU)); + volScalarField magS = mag(symm(gradU)); + forAll(inputFields[idxI], cellI) + { + inputFields[idxI][cellI] = magS[cellI] / magOmega[cellI]; + } + } + else if (inputName == "PoD") + { + daModel_.getTurbProdOverDestruct(inputFields[idxI]); + } + else if (inputName == "chiSA") + { +#ifndef SolidDASolver + const volScalarField& nuTilda = mesh_.thisDb().lookupObject("nuTilda"); + volScalarField nu = daModel_.getDATurbulenceModel().nu(); + forAll(inputFields[idxI], cellI) + { + inputFields[idxI][cellI] = nuTilda[cellI] / nu[cellI]; + } +#endif + } + else + { + FatalErrorIn("") << "inputName: " << inputName << " not supported. Options are: SoQ, PoD" << abort(FatalError); + } + } + + List> layerVals; + layerVals.setSize(nHiddenLayers); + for (label layerI = 0; layerI < nHiddenLayers; layerI++) + { + label nNeurons = hiddenLayerNeurons_[layerI]; + layerVals[layerI].setSize(nNeurons); + } + + forAll(mesh_.cells(), cellI) + { + label counterI = 0; + + for (label layerI = 0; layerI < nHiddenLayers; layerI++) + { + label nNeurons = hiddenLayerNeurons_[layerI]; + forAll(layerVals[layerI], neuronI) + { + layerVals[layerI][neuronI] = 0.0; + } + for (label neuronI = 0; neuronI < nNeurons; neuronI++) + { + if (layerI == 0) + { + forAll(inputNames_, neuronJ) + { + // weighted sum + layerVals[layerI][neuronI] += inputFields[neuronJ][cellI] * parameters_[counterI]; + counterI++; + } + } + else + { + forAll(layerVals[layerI - 1], neuronJ) + { + // weighted sum + layerVals[layerI][neuronI] += layerVals[layerI - 1][neuronJ] * parameters_[counterI]; + counterI++; + } + } + // bias + layerVals[layerI][neuronI] += parameters_[counterI]; + counterI++; + // activation function + layerVals[layerI][neuronI] = 1 / (1 + exp(-layerVals[layerI][neuronI])); + } + } + // final output layer, we have only one output + scalar outputVal = 0.0; + forAll(layerVals[nHiddenLayers - 1], neuronJ) + { + // weighted sum + outputVal += layerVals[nHiddenLayers - 1][neuronJ] * parameters_[counterI]; + counterI++; + } + // bias + outputVal += parameters_[counterI]; + + outputField[cellI] = outputScale_ * (outputVal + outputShift_); + } + } + else + { + FatalErrorIn("") << "modelType_: " << modelType_ << " not supported. Options are: neuralNetwork" << abort(FatalError); + } +} + +label DARegression::nParameters() +{ + /* + Description: + get the number of parameters + */ + + if (modelType_ == "neuralNetwork") + { + label nHiddenLayers = hiddenLayerNeurons_.size(); + label nInputs = inputNames_.size(); + + // add weights + // input + label nParameters = nInputs * hiddenLayerNeurons_[0]; + // hidden layers + for (label layerI = 1; layerI < hiddenLayerNeurons_.size(); layerI++) + { + nParameters += hiddenLayerNeurons_[layerI] * hiddenLayerNeurons_[layerI - 1]; + } + // output + nParameters += hiddenLayerNeurons_[nHiddenLayers - 1] * 1; + + // add biases + // add hidden layers + for (label layerI = 0; layerI < hiddenLayerNeurons_.size(); layerI++) + { + nParameters += hiddenLayerNeurons_[layerI]; + } + // add output layer + nParameters += 1; + + return nParameters; + } +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/adjoint/DARegression/DARegression.H b/src/adjoint/DARegression/DARegression.H new file mode 100755 index 00000000..dbe8df4b --- /dev/null +++ b/src/adjoint/DARegression/DARegression.H @@ -0,0 +1,110 @@ +/*---------------------------------------------------------------------------*\ + + DAFoam : Discrete Adjoint with OpenFOAM + Version : v3 + + Description: + Regression model + +\*---------------------------------------------------------------------------*/ + +#ifndef DARegression_H +#define DARegression_H + +#include "fvOptions.H" +#include "surfaceFields.H" +#include "DAOption.H" +#include "DAUtility.H" +#include "DAModel.H" +#include "globalIndex.H" +#include "DAMacroFunctions.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class DARegression Declaration +\*---------------------------------------------------------------------------*/ + +class DARegression +{ + +private: + /// Disallow default bitwise copy construct + DARegression(const DARegression&); + + /// Disallow default bitwise assignment + void operator=(const DARegression&); + +protected: + /// Foam::fvMesh object + const fvMesh& mesh_; + + /// Foam::DAOption object + const DAOption& daOption_; + + /// DAModel object + const DAModel& daModel_; + + /// the type of regression model + word modelType_; + + /// a list of words for the inputs + wordList inputNames_; + + /// a list of words for the outputs + word outputName_; + + /// number of neurons hidden layers of the neural network + labelList hiddenLayerNeurons_; + + /// we can shift the output. we always shift before scaling it. + scalar outputShift_; + + /// we can scale the output. we always shift before scaling it. + scalar outputScale_; + + /// the parameters for the regression model + scalarList parameters_; + +public: + /// Constructors + DARegression( + const fvMesh& mesh, + const DAOption& daOption, + const DAModel& daModel); + + /// Destructor + virtual ~DARegression() + { + } + + // Members + + void compute(); + + label nParameters(); + + scalar getParameter(label idxI) + { + return parameters_[idxI]; + } + + void setParameter(label idxI, scalar val) + { + parameters_[idxI] = val; + } + +}; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C b/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C index 99331025..dd4db875 100755 --- a/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C +++ b/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.C @@ -99,6 +99,8 @@ void DAPimpleFoam::initSolver() break; } } + + hasRegModel_ = daOptionPtr_->getAllOptions().subDict("regressionModel").getLabel("active"); } label DAPimpleFoam::solvePrimal( @@ -175,6 +177,16 @@ label DAPimpleFoam::solvePrimal( scalar deltaT = runTime.deltaT().value(); label nInstances = round(endTime / deltaT); + // check if the parameters are set in the Python layer + if (hasRegModel_) + { + scalar testVal = daRegressionPtr_->getParameter(0); + if (testVal > 1e15) + { + FatalErrorIn("") << "regressionModel is active but the parameter values are not set!" << abort(FatalError); + } + } + // main loop for (label iter = 1; iter <= nInstances; iter++) { @@ -208,6 +220,12 @@ label DAPimpleFoam::solvePrimal( #include "pEqnPimple.H" } + if (hasRegModel_) + { + // if the regression model is active, update the output field value at each iteration + daRegressionPtr_->compute(); + } + laminarTransport.correct(); daTurbulenceModelPtr_->correct(pimplePrintToScreen); } diff --git a/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.H b/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.H index 170d634b..129ea3ef 100755 --- a/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.H +++ b/src/adjoint/DASolver/DAPimpleFoam/DAPimpleFoam.H @@ -101,6 +101,9 @@ protected: /// whether to write mesh for the reduceIO label reduceIOWriteMesh_ = 0; + /// whether we have a regression model to update the field + label hasRegModel_ = 0; + public: TypeName("DAPimpleFoam"); // Constructors diff --git a/src/adjoint/DASolver/DASolver.C b/src/adjoint/DASolver/DASolver.C index 57a9290c..5853a2cc 100644 --- a/src/adjoint/DASolver/DASolver.C +++ b/src/adjoint/DASolver/DASolver.C @@ -43,7 +43,8 @@ DASolver::DASolver( daFieldPtr_(nullptr), daCheckMeshPtr_(nullptr), daLinearEqnPtr_(nullptr), - daResidualPtr_(nullptr) + daResidualPtr_(nullptr), + daRegressionPtr_(nullptr) #ifdef CODI_AD_REVERSE , globalADTape_(codi::RealReverse::getTape()) @@ -5866,6 +5867,125 @@ void DASolver::calcdRdThermalTPsiAD( #endif } +void DASolver::calcdRdRegParTPsiAD( + const double* volCoords, + const double* states, + const double* parameters, + const double* seeds, + double* product) +{ +#ifdef CODI_AD_REVERSE + /* + Description: + Compute the matrix-vector products [dR/dRegParameters]^T*Psi using reverse-mode AD + + Input: + + xvVec: the volume mesh coordinate vector + + wVec: the state variable vector + + parameters: the parameters for the regression model + + psi: the vector to multiply dRdXv + + Output: + prodVec: the matrix-vector products [dR/dRegParameters]^T*Psi + */ + + Info << "Calculating [dRdRegPar]^T * Psi 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]; + } + + this->updateOFMesh(volCoordsArray); + this->updateOFField(statesArray); + + // update the OpenFOAM variables and reset their seeds (gradient part) to zeros + this->resetOFSeeds(); + // reset the AD tape + this->globalADTape_.reset(); + // 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]); + } + + // calculate outputs + this->regressionModelCompute(); + + // compute residuals + this->updateStateBoundaryConditions(); + this->calcResiduals(); + + // register outputs + this->registerResidualOutput4AD(); + + // stop recording + this->globalADTape_.setPassive(); + + // set seeds to the outputs + this->assignSeeds2ResidualGradient(seeds); + + // now calculate the reverse matrix-vector product + this->globalADTape_.evaluate(); + + // get the matrix-vector product from the inputs + for (label i = 0; i < nParameters; i++) + { + product[i] = parametersArray[i].getGradient(); + } + + // clean up AD + 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->regressionModelCompute(); + this->updateStateBoundaryConditions(); + this->calcResiduals(); + + delete[] parametersArray; + delete[] volCoordsArray; + delete[] statesArray; + +#endif +} + void DASolver::calcdRdXvTPsiAD( const Vec xvVec, const Vec wVec, diff --git a/src/adjoint/DASolver/DASolver.H b/src/adjoint/DASolver/DASolver.H index 56d5f8c4..7bead9ab 100644 --- a/src/adjoint/DASolver/DASolver.H +++ b/src/adjoint/DASolver/DASolver.H @@ -33,6 +33,7 @@ #include "DAField.H" #include "DAPartDeriv.H" #include "DALinearEqn.H" +#include "DARegression.H" #include "volPointInterpolation.H" #include "IOMRFZoneListDF.H" #include "interpolateSplineXY.H" @@ -99,6 +100,9 @@ protected: /// DAStateInfo pointer autoPtr daStateInfoPtr_; + /// DARegression pointer + autoPtr daRegressionPtr_; + /// the stateInfo_ list from DAStateInfo object HashTable stateInfo_; @@ -1127,6 +1131,38 @@ public: return forwardADDerivVal_[objFuncName]; } + /// get the number of regression model parameters + label getNRegressionParameters() + { + return daRegressionPtr_->nParameters(); + } + + /// get the regression parameter + scalar getRegressionParameter(const label idxI) + { + return daRegressionPtr_->getParameter(idxI); + } + + /// set the regression parameter + void setRegressionParameter(const label idxI, scalar val) + { + daRegressionPtr_->setParameter(idxI, val); + } + + /// call the compute method of the regression model + void regressionModelCompute() + { + daRegressionPtr_->compute(); + } + + /// calculte [dR/dRegParameters]^T * psi + void calcdRdRegParTPsiAD( + const double* volCoords, + const double* states, + const double* parameters, + const double* seeds, + double* product); + /// update the primal state boundary condition based on the primalBC dict void setPrimalBoundaryConditions(const label printInfo = 1); diff --git a/src/adjoint/Make/files_Incompressible b/src/adjoint/Make/files_Incompressible index aec2dcf2..072f4eb8 100755 --- a/src/adjoint/Make/files_Incompressible +++ b/src/adjoint/Make/files_Incompressible @@ -27,6 +27,8 @@ DAModel/DARadiationModel/DAP1.C DAModel/DAModel.C +DARegression/DARegression.C + DAStateInfo/DAStateInfo.C DAStateInfo/DAStateInfoSimpleFoam.C DAStateInfo/DAStateInfoPisoFoam.C diff --git a/src/include/createAdjointIncompressible.H b/src/include/createAdjointIncompressible.H index 065ffb6e..92646b2b 100755 --- a/src/include/createAdjointIncompressible.H +++ b/src/include/createAdjointIncompressible.H @@ -27,6 +27,13 @@ daTurbulenceModelPtr_.reset(DATurbulenceModel::New(turbModelName, mesh, daOption daModelPtr_.reset(new DAModel(mesh, daOptionPtr_())); +label hasRegModel = daOptionPtr_->getAllOptions().subDict("regressionModel").getLabel("active"); +if (hasRegModel) +{ + Info << "Init DARegression " << endl; + daRegressionPtr_.reset(new DARegression(mesh, daOptionPtr_(), daModelPtr_())); +} + word solverName = daOptionPtr_->getOption("solverName"); daStateInfoPtr_.reset(DAStateInfo::New(solverName, mesh, daOptionPtr_(), daModelPtr_())); stateInfo_ = daStateInfoPtr_->getStateInfo(); diff --git a/src/include/setForwardADSeeds.H b/src/include/setForwardADSeeds.H index fd993a0e..7228ad4e 100755 --- a/src/include/setForwardADSeeds.H +++ b/src/include/setForwardADSeeds.H @@ -371,6 +371,19 @@ if (daOptionPtr_->getAllOptions().subDict("useAD").getWord("mode") == "forward") } } } + else if (type == "RegPar") + { + // for ACTD derivs + scalar seed = 0.0; + + label seedIndex = daOptionPtr_->getAllOptions().subDict("useAD").getLabel("seedIndex"); + + seed = daRegressionPtr_->getParameter(seedIndex); + + seed.setGradient(1.0); + + daRegressionPtr_->setParameter(seedIndex, seed); + } } #endif diff --git a/src/pyDASolvers/DASolvers.H b/src/pyDASolvers/DASolvers.H index 55694ff6..bf280fce 100755 --- a/src/pyDASolvers/DASolvers.H +++ b/src/pyDASolvers/DASolvers.H @@ -582,6 +582,36 @@ public: << abort(FatalError); } + /// get the number of regression model parameters + label getNRegressionParameters() + { + return DASolverPtr_->getNRegressionParameters(); + } + + /// call the compute method of the regression model + void regressionModelCompute() + { + DASolverPtr_->regressionModelCompute(); + } + + /// set the regression parameter + void setRegressionParameter(label idxI, double val) + { + scalar val1 = val; + DASolverPtr_->setRegressionParameter(idxI, val1); + } + + /// calculte [dR/dRegParameters]^T * psi + void calcdRdRegParTPsiAD( + const double* volCoords, + const double* states, + const double* parameters, + const double* seeds, + double* product) + { + DASolverPtr_->calcdRdRegParTPsiAD(volCoords, states, parameters, seeds, product); + } + void calcdRdThermalTPsiAD( const double* volCoords, const double* states, diff --git a/src/pyDASolvers/pyDASolvers.pyx b/src/pyDASolvers/pyDASolvers.pyx index 79916961..055e1933 100755 --- a/src/pyDASolvers/pyDASolvers.pyx +++ b/src/pyDASolvers/pyDASolvers.pyx @@ -75,6 +75,7 @@ cdef extern from "DASolvers.H" namespace "Foam": void calcdRdFieldTPsiAD(PetscVec, PetscVec, PetscVec, char *, PetscVec) void calcdFdFieldAD(PetscVec, PetscVec, char *, char *, PetscVec) void calcdRdThermalTPsiAD(double *, double *, double *, double *, double *) + void calcdRdRegParTPsiAD(double *, double *, double *, double *, double *) void calcdRdWOldTPsiAD(int, PetscVec, PetscVec) void convertMPIVec2SeqVec(PetscVec, PetscVec) void syncDAOptionToActuatorDVs() @@ -103,6 +104,9 @@ cdef extern from "DASolvers.H" namespace "Foam": void getThermal(double *, double *, double *) void getThermalAD(char *, double *, double *, double *, double *) void setThermal(double *) + int getNRegressionParameters() + void setRegressionParameter(int, double) + void regressionModelCompute() void getOFField(char *, char *, PetscVec) void getAcousticData(PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, PetscVec, char*) void printAllOptions() @@ -471,6 +475,41 @@ cdef class pyDASolvers: assert len(thermal) == self.getNCouplingFaces() * 2, "invalid array size!" cdef double *thermal_data = thermal.data self._thisptr.setThermal(thermal_data) + + def getNRegressionParameters(self): + return self._thisptr.getNRegressionParameters() + + def setRegressionParameter(self, idx, val): + self._thisptr.setRegressionParameter(idx, val) + + def regressionModelCompute(self): + self._thisptr.regressionModelCompute() + + def calcdRdRegParTPsiAD(self, + np.ndarray[double, ndim=1, mode="c"] volCoords, + np.ndarray[double, ndim=1, mode="c"] states, + np.ndarray[double, ndim=1, mode="c"] parameters, + np.ndarray[double, ndim=1, mode="c"] seeds, + np.ndarray[double, ndim=1, mode="c"] product): + + assert len(volCoords) == self.getNLocalPoints() * 3, "invalid array size!" + assert len(states) == self.getNLocalAdjointStates(), "invalid array size!" + assert len(parameters) == self.getNRegressionParameters(), "invalid array size!" + assert len(seeds) == self.getNLocalAdjointStates(), "invalid array size!" + assert len(product) == self.getNRegressionParameters(), "invalid array size!" + + cdef double *volCoords_data = volCoords.data + cdef double *states_data = states.data + cdef double *parameters_data = parameters.data + cdef double *seeds_data = seeds.data + cdef double *product_data = product.data + + self._thisptr.calcdRdRegParTPsiAD( + volCoords_data, + states_data, + parameters_data, + seeds_data, + product_data) def getAcousticData(self, Vec x, Vec y, Vec z, Vec nX, Vec nY, Vec nZ, Vec a, Vec fX, Vec fY, Vec fZ, groupName): self._thisptr.getAcousticData(x.vec, y.vec, z.vec, nX.vec, nY.vec, nZ.vec, a.vec, fX.vec, fY.vec, fZ.vec, groupName) From e4197e66479b16b5eff44c4a8190e8eee03cb8dd Mon Sep 17 00:00:00 2001 From: Ping He Date: Thu, 28 Dec 2023 10:31:12 -0600 Subject: [PATCH 2/8] Called correctBC for regModel. --- src/adjoint/DARegression/DARegression.C | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/adjoint/DARegression/DARegression.C b/src/adjoint/DARegression/DARegression.C index f77c8391..5df08a13 100755 --- a/src/adjoint/DARegression/DARegression.C +++ b/src/adjoint/DARegression/DARegression.C @@ -171,6 +171,8 @@ void DARegression::compute() outputField[cellI] = outputScale_ * (outputVal + outputShift_); } + + outputField.correctBoundaryConditions(); } else { From 41ab432e00e19efe3395cd66f333165906e34102 Mon Sep 17 00:00:00 2001 From: Ping He Date: Thu, 28 Dec 2023 10:31:46 -0600 Subject: [PATCH 3/8] Fixed a missing init bug. --- src/adjoint/DASolver/DASolver.H | 2 ++ src/adjoint/Make/files_Compressible | 2 ++ src/include/createAdjointCompressible.H | 7 +++++++ 3 files changed, 11 insertions(+) diff --git a/src/adjoint/DASolver/DASolver.H b/src/adjoint/DASolver/DASolver.H index 7bead9ab..e17ed8f7 100644 --- a/src/adjoint/DASolver/DASolver.H +++ b/src/adjoint/DASolver/DASolver.H @@ -1131,6 +1131,7 @@ public: return forwardADDerivVal_[objFuncName]; } +#ifndef SolidDASolver /// get the number of regression model parameters label getNRegressionParameters() { @@ -1162,6 +1163,7 @@ public: const double* parameters, const double* seeds, double* product); +#endif /// update the primal state boundary condition based on the primalBC dict void setPrimalBoundaryConditions(const label printInfo = 1); diff --git a/src/adjoint/Make/files_Compressible b/src/adjoint/Make/files_Compressible index 1901dcca..bdd91278 100755 --- a/src/adjoint/Make/files_Compressible +++ b/src/adjoint/Make/files_Compressible @@ -26,6 +26,8 @@ DAModel/DARadiationModel/DAP1.C DAModel/DAModel.C +DARegression/DARegression.C + DAStateInfo/DAStateInfo.C DAStateInfo/DAStateInfoRhoSimpleFoam.C DAStateInfo/DAStateInfoRhoSimpleCFoam.C diff --git a/src/include/createAdjointCompressible.H b/src/include/createAdjointCompressible.H index f42a4267..a1363121 100755 --- a/src/include/createAdjointCompressible.H +++ b/src/include/createAdjointCompressible.H @@ -27,6 +27,13 @@ daTurbulenceModelPtr_.reset(DATurbulenceModel::New(turbModelName, mesh, daOption daModelPtr_.reset(new DAModel(mesh, daOptionPtr_())); +label hasRegModel = daOptionPtr_->getAllOptions().subDict("regressionModel").getLabel("active"); +if (hasRegModel) +{ + Info << "Init DARegression " << endl; + daRegressionPtr_.reset(new DARegression(mesh, daOptionPtr_(), daModelPtr_())); +} + word solverName = daOptionPtr_->getOption("solverName"); daStateInfoPtr_.reset(DAStateInfo::New(solverName, mesh, daOptionPtr_(), daModelPtr_())); stateInfo_ = daStateInfoPtr_->getStateInfo(); From af630c335c9b432580f00b04f21a207444abc7cd Mon Sep 17 00:00:00 2001 From: Ping He Date: Thu, 28 Dec 2023 11:39:10 -0600 Subject: [PATCH 4/8] Added an option to not time vol for the volSum obj. --- src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C | 14 ++++++++++++-- src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.H | 3 +++ 2 files changed, 15 insertions(+), 2 deletions(-) diff --git a/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C b/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C index 0ec19b99..3c9b3718 100755 --- a/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C +++ b/src/adjoint/DAObjFunc/DAObjFuncVariableVolSum.C @@ -51,6 +51,8 @@ DAObjFuncVariableVolSum::DAObjFuncVariableVolSum( objFuncDict_.readEntry