DASolver.H

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
 
    DAFoam  : Discrete Adjoint with OpenFOAM
    Version : v3
 
    Description:
        Solver class that calls primal and adjoint solvers,
        and compute the total derivatives
 
\*---------------------------------------------------------------------------*/
 
#ifndef DASolver_H
#define DASolver_H
 
#include <petscksp.h>
#include "Python.h"
#include "fvCFD.H"
#include "fvMesh.H"
#include "runTimeSelectionTables.H"
#include "OFstream.H"
#include "functionObjectList.H"
#include "fvOptions.H"
#include "DAUtility.H"
#include "DACheckMesh.H"
#include "DAOption.H"
#include "DAStateInfo.H"
#include "DAModel.H"
#include "DAIndex.H"
#include "DAObjFunc.H"
#include "DAJacCon.H"
#include "DAColoring.H"
#include "DAResidual.H"
#include "DAField.H"
#include "DAPartDeriv.H"
#include "DALinearEqn.H"
#include "volPointInterpolation.H"
#include "IOMRFZoneListDF.H"
#include "interpolateSplineXY.H"
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
namespace Foam
{
 
/*---------------------------------------------------------------------------*\
                    Class DASolver Declaration
\*---------------------------------------------------------------------------*/
 
class DASolver
{
 
private:
    DASolver(const DASolver&);
 
    void operator=(const DASolver&);
 
protected:
    char* argsAll_;
 
    PyObject* pyOptions_;
 
    autoPtr<argList> argsPtr_;
 
    autoPtr<Time> runTimePtr_;
 
    autoPtr<fvMesh> meshPtr_;
 
    autoPtr<DAOption> daOptionPtr_;
 
    autoPtr<DAModel> daModelPtr_;
 
    autoPtr<DAIndex> daIndexPtr_;
 
    autoPtr<DAField> daFieldPtr_;
 
    UPtrList<DAObjFunc> daObjFuncPtrList_;
 
    autoPtr<DACheckMesh> daCheckMeshPtr_;
 
    autoPtr<DALinearEqn> daLinearEqnPtr_;
 
    autoPtr<DAResidual> daResidualPtr_;
 
    autoPtr<DAStateInfo> daStateInfoPtr_;
 
    HashTable<wordList> stateInfo_;
 
    wordList objFuncNames4Adj_;
 
    dictionary psiVecDict_;
 
    dictionary totalDerivDict_;
 
    autoPtr<OFstream> objFuncAvgHistFilePtr_;
 
    label nItersObjFuncAvg_ = -9999;
 
    scalarList avgObjFuncValues_;
 
    Mat dRdWTPC_;
 
    label nSolveAdjointCalls_ = 0;
 
    Mat dXvdFFDMat_;
 
    Vec FFD2XvSeedVec_;
 
    HashTable<PetscScalar> forwardADDerivVal_;
 
    scalar primalMinRes_ = 1e10;
 
    scalar prevPrimalSolTime_ = -1e10;
 
    template<class classType>
    void primalResidualControl(
        const SolverPerformance<classType>& solverP,
        const label printToScreen,
        const label printInterval,
        const word varName);
 
    label isPrintTime(
        const Time& runTime,
        const label printInterval) const;
 
    label checkResidualTol();
 
    void reduceStateResConLevel(
        const dictionary& maxResConLv4JacPCMat,
        HashTable<List<List<word>>>& stateResConInfo) const;
 
    void writeAssociatedFields();
 
    Mat dRdWTMF_;
 
    label globalADTape4dRdWTInitialized = 0;
 
    List<List<scalar>> stateAllInstances_;
 
    List<List<scalar>> stateBoundaryAllInstances_;
 
    List<dictionary> objFuncsAllInstances_;
 
    List<scalar> runTimeAllInstances_;
 
    List<label> runTimeIndexAllInstances_;
 
    label nTimeInstances_ = -9999;
 
    scalar periodicity_ = 0.0;
 
    scalar primalMinResTol_ = 0.0;
 
    label primalMinIters_ = -1;
 
    void saveTimeInstanceFieldHybrid(label& timeInstanceI);
 
    void saveTimeInstanceFieldTimeAccurate(label& timeInstanceI);
 
public:
    TypeName("DASolver");
 
    // Declare run-time constructor selection table
    declareRunTimeSelectionTable(
        autoPtr,
        DASolver,
        dictionary,
        (
            char* argsAll,
            PyObject* pyOptions),
        (argsAll, pyOptions));
 
    // Constructors
 
    //- Construct from components
    DASolver(
        char* argsAll,
        PyObject* pyOptions);
 
    // Selectors
 
    //- Return a reference to the selected model
    static autoPtr<DASolver> New(
        char* argsAll,
        PyObject* pyOptions);
 
    //- Destructor
    virtual ~DASolver()
    {
    }
 
    // Member functions
 
    virtual void initSolver() = 0;
 
    virtual label solvePrimal(
        const Vec xvVec,
        Vec wVec) = 0;
 
    virtual label runFPAdj(
        const Vec xvVec,
        const Vec wVec,
        Vec dFdW,
        Vec psi);
 
    void setTimeInstanceField(const label instanceI);
 
    scalar getTimeInstanceObjFunc(
        const label instanceI,
        const word objFuncName);
 
    void setTimeInstanceVar(
        const word mode,
        Mat stateMat,
        Mat stateBCMat,
        Vec timeVec,
        Vec timeIdxVec);
 
    void initOldTimes();
 
    void calcdRdWT(
        const Vec xvVec,
        const Vec wVec,
        const label isPC,
        Mat dRdWT);
 
    void calcdRdWTPsiAD(
        const Vec xvVec,
        const Vec wVec,
        const Vec psi,
        Vec dRdWTPsi);
 
    void calcdRdWTPsiAD(
        const label isInit,
        const Vec psi,
        Vec dRdWTPsi);
 
    void calcdFdW(
        const Vec xvVec,
        const Vec wVec,
        const word objFuncName,
        Vec dFdW);
 
    void createMLRKSP(
        const Mat jacMat,
        const Mat jacPCMat,
        KSP ksp);
 
    label solveLinearEqn(
        const KSP ksp,
        const Vec rhsVec,
        Vec solVec);
 
    void convertMPIVec2SeqVec(
        const Vec mpiVec,
        Vec seqVec);
 
    void updateOFField(const Vec wVec);
 
    void updateOFMesh(const Vec xvVec);
 
    void updateOFField(const scalar* volCoords);
 
    void updateOFMesh(const scalar* states);
 
    void resetOFSeeds();
 
    void calcdRdBC(
        const Vec xvVec,
        const Vec wVec,
        const word designVarName,
        Mat dRdBC);
 
    void calcdFdBC(
        const Vec xvVec,
        const Vec wVec,
        const word objFuncName,
        const word designVarName,
        Vec dFdBC);
 
    void calcdRdBCTPsiAD(
        const Vec xvVec,
        const Vec wVec,
        const Vec psi,
        const word designVarName,
        Vec dRdBCTPsi);
 
    void calcdFdBCAD(
        const Vec xvVec,
        const Vec wVec,
        const word objFuncName,
        const word designVarName,
        Vec dFdBC);
 
    void calcdRdAOA(
        const Vec xvVec,
        const Vec wVec,
        const word designVarName,
        Mat dRdAOA);
 
    void calcdFdAOA(
        const Vec xvVec,
        const Vec wVec,
        const word objFuncName,
        const word designVarName,
        Vec dFdAOA);
 
    void calcdRdAOATPsiAD(
        const Vec xvVec,
        const Vec wVec,
        const Vec psi,
        const word designVarName,
        Vec dRdAOATPsi);
 
    void calcdRdFFD(
        const Vec xvVec,
        const Vec wVec,
        const word designVarName,
        Mat dRdFFD);
 
    void calcdFdFFD(
        const Vec xvVec,
        const Vec wVec,
        const word objFuncName,
        const word designVarName,
        Vec dFdFFD);
 
    void calcdRdACT(
        const Vec xvVec,
        const Vec wVec,
        const word designVarName,
        const word designVarType,
        Mat dRdACT);
 
    void calcdFdACT(
        const Vec xvVec,
        const Vec wVec,
        const word objFuncName,
        const word designVarName,
        const word designVarType,
        Vec dFdACT);
 
    void calcdRdFieldTPsiAD(
        const Vec xvVec,
        const Vec wVec,
        const Vec psi,
        const word designVarName,
        Vec dRdFieldTPsi);
 
    void calcdFdFieldAD(
        const Vec xvVec,
        const Vec wVec,
        const word objFuncName,
        const word designVarName,
        Vec dFdField);
 
    void createMLRKSPMatrixFree(
        const Mat jacPCMat,
        KSP ksp);
 
    void calcdFdWAD(
        const Vec xvVec,
        const Vec wVec,
        const word objFuncName,
        Vec dFdW);
 
    void calcdRdThermalTPsiAD(
        const double* volCoords,
        const double* states,
        const double* thermal,
        const double* seeds,
        double* product);
 
    void calcdRdXvTPsiAD(
        const Vec xvVec,
        const Vec wVec,
        const Vec psi,
        Vec dRdXvTPsi);
 
    void calcdForcedXvAD(
        const Vec xvVec,
        const Vec wVec,
        const Vec fBarVec,
        Vec dForcedXv);
 
    void calcdAcousticsdXvAD(
        const Vec xvVec,
        const Vec wVec,
        const Vec fBarVec,
        Vec dForcedXv,
        word varName,
        word groupName);
 
    void calcdFdXvAD(
        const Vec xvVec,
        const Vec wVec,
        const word objFuncName,
        const word designVarName,
        Vec dFdXv);
 
    void calcdRdActTPsiAD(
        const Vec xvVec,
        const Vec wVec,
        const Vec psi,
        const word designVarName,
        Vec dRdActTPsi);
 
    void calcdForcedWAD(
        const Vec xvVec,
        const Vec wVec,
        const Vec fBarVec,
        Vec dForcedW);
 
    void calcdAcousticsdWAD(
        const Vec xvVec,
        const Vec wVec,
        const Vec fBarVec,
        Vec dForcedW,
        word varName,
        word groupName);
 
    void calcdFdACTAD(
        const Vec xvVec,
        const Vec wVec,
        const word objFuncName,
        const word designVarName,
        Vec dFdACT);
 
    void calcdRdWOldTPsiAD(
        const label oldTimeLevel,
        const Vec psi,
        Vec dRdWOldTPsi);
 
    void getCouplingPatchList(wordList& patchList, word groupName = "NONE");
 
    static PetscErrorCode dRdWTMatVecMultFunction(
        Mat dRdWT,
        Vec vecX,
        Vec vecY);
 
    void initializedRdWTMatrixFree(
        const Vec xvVec,
        const Vec wVec);
 
    void destroydRdWTMatrixFree();
 
    void registerStateVariableInput4AD(const label oldTimeLevel = 0);
 
    void registerFieldVariableInput4AD(
        const word fieldName,
        const word fieldType);
 
    void registerResidualOutput4AD();
 
    void registerForceOutput4AD(List<scalar>& fX, List<scalar>& fY, List<scalar>& fZ);
 
    void registerAcousticOutput4AD(List<scalar>& a);
 
    void assignVec2ResidualGradient(Vec vecX);
 
    void assignSeeds2ResidualGradient(const double* seeds);
 
    void assignVec2ForceGradient(Vec fBarVec, List<scalar>& fX, List<scalar>& fY, List<scalar>& fZ);
 
    void assignVec2AcousticGradient(Vec fBarVec, List<scalar>& a, label offset, label step);
 
    void assignStateGradient2Vec(
        Vec vecY,
        const label oldTimeLevel = 0);
 
    void assignFieldGradient2Vec(
        const word fieldName,
        const word fieldType,
        Vec vecY);
 
    void normalizeGradientVec(Vec vecY);
 
    void initializeGlobalADTape4dRdWT();
 
    label loop(Time& runTime);
 
    label getGlobalXvIndex(
        const label idxPoint,
        const label idxCoord) const
    {
        return daIndexPtr_->getGlobalXvIndex(idxPoint, idxCoord);
    }
 
    void ofField2StateVec(Vec stateVec) const
    {
        daFieldPtr_->ofField2StateVec(stateVec);
    }
 
    void stateVec2OFField(const Vec stateVec) const
    {
        daFieldPtr_->stateVec2OFField(stateVec);
    }
 
    void pointVec2OFMesh(const Vec xvVec) const
    {
        daFieldPtr_->pointVec2OFMesh(xvVec);
    }
 
    void ofMesh2PointVec(Vec xvVec) const
    {
        daFieldPtr_->ofMesh2PointVec(xvVec);
    }
 
    void resVec2OFResField(const Vec resVec) const
    {
        daFieldPtr_->resVec2OFResField(resVec);
    }
 
    void ofResField2ResVec(Vec resVec) const
    {
        daFieldPtr_->resVec2OFResField(resVec);
    }
 
    void writeMatrixBinary(
        const Mat matIn,
        const word prefix)
    {
        DAUtility::writeMatrixBinary(matIn, prefix);
    }
 
    void writeMatrixASCII(
        const Mat matIn,
        const word prefix)
    {
        DAUtility::writeMatrixASCII(matIn, prefix);
    }
 
    void readMatrixBinary(
        Mat matIn,
        const word prefix)
    {
        DAUtility::readMatrixBinary(matIn, prefix);
    }
 
    void writeVectorASCII(
        const Vec vecIn,
        const word prefix)
    {
        DAUtility::writeVectorASCII(vecIn, prefix);
    }
 
    void readVectorBinary(
        Vec vecIn,
        const word prefix)
    {
        DAUtility::readVectorBinary(vecIn, prefix);
    }
 
    void writeVectorBinary(
        const Vec vecIn,
        const word prefix)
    {
        DAUtility::writeVectorBinary(vecIn, prefix);
    }
 
    label getNLocalAdjointStates() const
    {
        return daIndexPtr_->nLocalAdjointStates;
    }
 
    label getNLocalAdjointBoundaryStates() const
    {
        return daIndexPtr_->nLocalAdjointBoundaryStates;
    }
 
    label getNLocalCells() const
    {
        return meshPtr_->nCells();
    }
 
    label getNLocalPoints() const
    {
        return meshPtr_->nPoints();
    }
 
    void setDAObjFuncList();
 
    void printAllObjFuncs();
 
    label checkMesh() const
    {
        return daCheckMeshPtr_->run();
    }
 
    scalar getObjFuncValue(const word objFuncName);
 
    void getForces(
        Vec fX,
        Vec fY,
        Vec fZ);
 
    void getAcousticData(
        Vec x,
        Vec y,
        Vec z,
        Vec nX,
        Vec nY,
        Vec nZ,
        Vec a,
        Vec fX,
        Vec fY,
        Vec fZ,
        word groupName);
 
    void getPatchInfo(
        label& nPoints,
        label& nFaces,
        List<word>& patchList);
 
    void getForcesInternal(
        List<scalar>& fX,
        List<scalar>& fY,
        List<scalar>& fZ,
        List<word>& patchList);
 
    label getNCouplingFaces();
 
    label getNCouplingPoints();
 
    void calcCouplingFaceCoords(
        const scalar* volCoords,
        scalar* surfCoords);
 
    void calcCouplingFaceCoordsAD(
        const double* volCoords,
        const double* seeds,
        double* product);
 
    void getThermal(
        const scalar* volCoords,
        const scalar* states,
        scalar* thermal);
 
    void getThermalAD(
        const word inputName,
        const double* volCoords,
        const double* states,
        const double* seeds,
        double* product);
 
    void setThermal(scalar* thermal);
 
    void getAcousticDataInternal(
        List<scalar>& x,
        List<scalar>& y,
        List<scalar>& z,
        List<scalar>& nX,
        List<scalar>& nY,
        List<scalar>& nZ,
        List<scalar>& a,
        List<scalar>& fX,
        List<scalar>& fY,
        List<scalar>& fZ,
        List<word>& patchList);
 
    void calcForceProfile(
        Vec center,
        Vec aForceL,
        Vec tForceL,
        Vec rDistL);
 
    void calcForceProfileInternal(
        fvMesh& mesh,
        const vector& center,
        scalarList& aForceL,
        scalarList& tForceL,
        scalarList& rDistL);
 
    void calcdForcedStateTPsiAD(
        const word mode,
        Vec xvVec,
        Vec stateVec,
        Vec psiVec,
        Vec prodVec);
 
    void calcFvSourceInternal(
        const word propName,
        const scalarField& aForce,
        const scalarField& tForce,
        const scalarField& rDistList,
        const scalarList& targetForce,
        const vector& center,
        volVectorField& fvSource);
 
    void calcFvSource(
        const word propName,
        Vec aForce,
        Vec tForce,
        Vec rDist,
        Vec targetForce,
        Vec center,
        Vec fvSource);
 
    void calcdFvSourcedInputsTPsiAD(
        const word propName,
        const word mode,
        Vec aForce,
        Vec tForce,
        Vec rDist,
        Vec targetForce,
        Vec center,
        Vec psi,
        Vec dFvSource);
 
    void printAllOptions()
    {
        Info << "DAFoam option dictionary: ";
        Info << daOptionPtr_->getAllOptions() << endl;
    }
 
    void calcPrimalResidualStatistics(
        const word mode,
        const label writeRes = 0);
 
    void setdXvdFFDMat(const Mat dXvdFFDMat);
 
    void setFFD2XvSeedVec(Vec vecIn);
 
    void updateDAOption(PyObject* pyOptions)
    {
        daOptionPtr_->updateDAOption(pyOptions);
    }
 
    scalar getPrevPrimalSolTime()
    {
        return prevPrimalSolTime_;
    }
 
    void setFieldValue4GlobalCellI(
        const word fieldName,
        const scalar val,
        const label globalCellI,
        const label compI = 0);
 
    void setFieldValue4LocalCellI(
        const word fieldName,
        const scalar val,
        const label localCellI,
        const label compI = 0);
 
    void updateBoundaryConditions(
        const word fieldName,
        const word fieldType);
 
    void syncDAOptionToActuatorDVs()
    {
        DAFvSource& fvSource = const_cast<DAFvSource&>(
            meshPtr_->thisDb().lookupObject<DAFvSource>("DAFvSource"));
        fvSource.syncDAOptionToActuatorDVs();
    }
 
    const fvMesh& getMesh()
    {
        return meshPtr_();
    }
 
    const Time& getRunTime()
    {
        return runTimePtr_();
    }
 
    const DAOption& getDAOption()
    {
        return daOptionPtr_();
    }
 
    const DAStateInfo& getDAStateInfo()
    {
        return daStateInfoPtr_();
    }
 
    const DAIndex& getDAIndex()
    {
        return daIndexPtr_();
    }
 
    const DAModel& getDAModel()
    {
        return daModelPtr_();
    }
 
    const DAResidual& getDAResidual()
    {
        return daResidualPtr_();
    }
 
    const DAField& getDAField()
    {
        return daFieldPtr_();
    }
 
    const DALinearEqn& getDALinearEqn()
    {
        return daLinearEqnPtr_();
    }
 
    const DACheckMesh& getDACheckMesh()
    {
        return daCheckMeshPtr_();
    }
 
    PetscScalar getForwardADDerivVal(const word objFuncName)
    {
        return forwardADDerivVal_[objFuncName];
    }
 
    void setPrimalBoundaryConditions(const label printInfo = 1);
 
    void calcResidualVec(Vec resVec);
 
    void writeFailedMesh();
 
    void initTensorFlowFuncs(
        pyComputeInterface computeInterface, 
        void* compute,
        pyJacVecProdInterface jacVecProdInterface,
        void* jacVecProd)
    {
        DAUtility::pyCalcBetaInterface = computeInterface;
        DAUtility::pyCalcBeta = compute;
 
        DAUtility::pyCalcBetaJacVecProdInterface = jacVecProdInterface;
        DAUtility::pyCalcBetaJacVecProd = jacVecProd;
    }
 
#ifdef CODI_AD_REVERSE
 
    codi::RealReverse::Tape& globalADTape_;
 
#endif
};
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
template<class classType>
void DASolver::primalResidualControl(
    const SolverPerformance<classType>& solverP,
    const label printToScreen,
    const label printInterval,
    const word varName)
{
    /*
    Description:
        Setup maximal residual control and print the residual as needed
    */
    if (mag(solverP.initialResidual()) < primalMinRes_)
    {
        primalMinRes_ = mag(solverP.initialResidual());
    }
    if (printToScreen)
    {
        Info << varName << " Initial residual: " << solverP.initialResidual() << endl
             << varName << "   Final residual: " << solverP.finalResidual() << endl;
    }
}
 
} // End namespace Foam
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
#endif
 
// ************************************************************************* //