DASolver.H

Go to the documentation of this file.

/*---------------------------------------------------------------------------*\
 
    DAFoam  : Discrete Adjoint with OpenFOAM
    Version : v4
 
    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 "DAFunction.H"
#include "DAJacCon.H"
#include "DAColoring.H"
#include "DAResidual.H"
#include "DAField.H"
#include "DAPartDeriv.H"
#include "DALinearEqn.H"
#include "DARegression.H"
#include "volPointInterpolation.H"
#include "IOMRFZoneListDF.H"
#include "interpolateSplineXY.H"
#include "DAInput.H"
#include "DAOutput.H"
#include "DAGlobalVar.H"
#include "DATimeOp.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<DAFunction> daFunctionPtrList_;
 
    UPtrList<DATimeOp> daTimeOpPtrList_;
 
    autoPtr<DACheckMesh> daCheckMeshPtr_;
 
    autoPtr<DALinearEqn> daLinearEqnPtr_;
 
    autoPtr<DAResidual> daResidualPtr_;
 
    autoPtr<DAStateInfo> daStateInfoPtr_;
 
    autoPtr<DARegression> daRegressionPtr_;
 
    autoPtr<DAGlobalVar> daGlobalVarPtr_;
 
    autoPtr<pointField> points0Ptr_;
 
    HashTable<wordList> stateInfo_;
 
    scalar prevPrimalSolTime_ = -1e10;
 
    label isPrintTime(
        const Time& runTime,
        const label printInterval) const;
 
    label checkPrimalFailure();
 
    void reduceStateResConLevel(
        const dictionary& maxResConLv4JacPCMat,
        HashTable<List<List<word>>>& stateResConInfo) const;
 
    void writeAssociatedFields();
 
    Mat dRdWTMF_;
 
    label globalADTape4dRdWTInitialized = 0;
 
    scalar primalMinResTol_ = 0.0;
 
    label printToScreen_ = 0;
 
    label regModelFail_ = 0;
 
    label printInterval_ = 100;
 
    label primalMinIters_ = -1;
 
    label printIntervalUnsteady_ = 1;
 
    List<scalarList> functionTimeSteps_;
 
    label primalFinalTimeIndex_;
 
    HashTable<scalar> initStateVals_;
 
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() = 0;
 
    virtual label runFPAdj(
        Vec dFdW,
        Vec psi);
 
    virtual label solveAdjointFP(
        Vec dFdW,
        Vec psi);
 
    void setTime(scalar time, label timeIndex)
    {
        runTimePtr_->setTime(time, timeIndex);
    }
 
    label getDdtSchemeOrder()
    {
        const fvSchemes& myFvSchemes = meshPtr_->thisDb().lookupObject<fvSchemes>("fvSchemes");
 
        word ddtSchemeName = myFvSchemes.subDict("ddtSchemes").getWord("default");
 
        if (ddtSchemeName == "steadyState")
        {
            return 0;
        }
        if (ddtSchemeName == "Euler")
        {
            return 1;
        }
        else if (ddtSchemeName == "backward")
        {
            return 2;
        }
        else
        {
            FatalErrorIn("") << "ddtScheme " << ddtSchemeName << " not supported! Options: steadyState, Euler, or backward"
                             << abort(FatalError);
            return -1;
        }
    }
 
    void printElapsedTime(const Time& runTime, const label printToScreen)
    {
        if (printToScreen)
        {
            Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
                 << "  ClockTime = " << runTime.elapsedClockTime() << " s"
                 << nl << endl;
        }
    }
 
    void calcdRdWT(
        const label isPC,
        Mat dRdWT);
 
    void updateKSPPCMat(
        Mat PCMat,
        KSP ksp);
 
    label solveLinearEqn(
        const KSP ksp,
        const Vec rhsVec,
        Vec solVec);
 
    void updateOFFields(const scalar* states);
 
    void updateOFMesh(const scalar* volCoords);
 
    void getOFFields(scalar* states)
    {
        daFieldPtr_->ofField2State(states);
    }
 
    void getOFField(
        const word fieldName,
        const word fieldType,
        double* field);
 
    void getOFMeshPoints(double* points);
 
    label getInputSize(
        const word inputName,
        const word inputType);
 
    label getOutputSize(
        const word outputName,
        const word outputType);
 
    label getInputDistributed(
        const word inputName,
        const word inputType);
 
    label getOutputDistributed(
        const word outputName,
        const word outputType);
 
    void calcOutput(
        const word outputName,
        const word outputType,
        double* output);
 
    void runColoring();
 
    void calcJacTVecProduct(
        const word inputName,
        const word inputType,
        const double* input,
        const word outputName,
        const word outputType,
        const double* seed,
        double* product);
 
    void setSolverInput(
        const word inputName,
        const word inputType,
        const int inputSize,
        const double* input,
        const double* seed);
 
    void createMLRKSPMatrixFree(
        const Mat jacPCMat,
        KSP ksp);
 
    void calcdRdWOldTPsiAD(
        const label oldTimeLevel,
        const double* psi,
        double* dRdWOldTPsi);
 
    void calcCouplingFaceCoords(
        const scalar* volCoords,
        scalar* surfCoords);
 
    void getCouplingPatchList(wordList& patchList, word groupName = "NONE");
 
    static PetscErrorCode dRdWTMatVecMultFunction(
        Mat dRdWT,
        Vec vecX,
        Vec vecY);
 
    void initializedRdWTMatrixFree();
 
    void destroydRdWTMatrixFree();
 
    void registerStateVariableInput4AD(const label oldTimeLevel = 0);
 
    void deactivateStateVariableInput4AD(const label oldTimeLevel = 0);
 
    void registerResidualOutput4AD();
 
    void assignVec2ResidualGradient(const double* vecX);
 
    void assignStateGradient2Vec(
        double* vecY,
        const label oldTimeLevel = 0);
 
    void normalizeGradientVec(double* vecY);
 
    void normalizeJacTVecProduct(
        const word inputName,
        double* product);
 
    void initializeGlobalADTape4dRdWT();
 
    label hasVolCoordInput();
 
    void initDynamicMesh();
 
    label loop(Time& runTime);
 
    void updateInputFieldUnsteady();
 
    void initInputFieldUnsteady();
 
    void meanStatesToStates();
 
    label getGlobalXvIndex(
        const label idxPoint,
        const label idxCoord) const
    {
        return daIndexPtr_->getGlobalXvIndex(idxPoint, idxCoord);
    }
 
    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 setDAFunctionList();
 
    void calcAllFunctions(label print = 0);
 
    double getTimeOpFuncVal(const word functionName);
 
    label getFunctionListIndex(const word functionName)
    {
        forAll(daFunctionPtrList_, idxI)
        {
            DAFunction& daFunction = daFunctionPtrList_[idxI];
            word functionName1 = daFunction.getFunctionName();
            if (functionName1 == functionName)
            {
                return idxI;
            }
        }
    }
 
    label checkMesh() const
    {
        return daCheckMeshPtr_->run();
    }
 
    scalar getdFScaling(const word functionName, const label timeIdx);
 
    void printAllOptions()
    {
        Info << "DAFoam option dictionary: ";
        Info << daOptionPtr_->getAllOptions() << endl;
    }
 
    void calcPrimalResidualStatistics(
        const word mode,
        const label writeRes = 0);
 
    void updateDAOption(PyObject* pyOptions)
    {
        daOptionPtr_->updateDAOption(pyOptions);
    }
 
    scalar getPrevPrimalSolTime()
    {
        return prevPrimalSolTime_;
    }
 
    void updateBoundaryConditions(
        const word fieldName,
        const word fieldType);
 
    void updateStateBoundaryConditions();
 
    void calcResiduals(label isPC = 0);
 
    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_();
    }
 
    label getNRegressionParameters(word modelName)
    {
        return daRegressionPtr_->nParameters(modelName);
    }
 
    scalar getRegressionParameter(word modelName, const label idxI)
    {
        return daRegressionPtr_->getParameter(modelName, idxI);
    }
 
    void setRegressionParameter(word modelName, const label idxI, scalar val)
    {
        daRegressionPtr_->setParameter(modelName, idxI, val);
    }
 
    void regressionModelCompute()
    {
        daRegressionPtr_->compute();
    }
 
    void setPrimalBoundaryConditions(const label printInfo = 1);
 
    void writeFailedMesh();
 
    void readStateVars(
        scalar timeVal,
        label oldTimeLevel = 0);
 
    void readMeshPoints(const scalar timeVal);
 
    void writeMeshPoints(const double* points, const scalar timeVal);
 
    void calcPCMatWithFvMatrix(Mat PCMat, const label turbOnly = 0);
 
    void initTensorFlowFuncs(
        pyComputeInterface computeInterface,
        void* compute,
        pyJacVecProdInterface jacVecProdInterface,
        void* jacVecProd,
        pySetCharInterface setModelNameInterface,
        void* setModelName)
    {
        DAUtility::pyCalcBetaInterface = computeInterface;
        DAUtility::pyCalcBeta = compute;
 
        DAUtility::pyCalcBetaJacVecProdInterface = jacVecProdInterface;
        DAUtility::pyCalcBetaJacVecProd = jacVecProd;
 
        DAUtility::pySetModelNameInterface = setModelNameInterface;
        DAUtility::pySetModelName = setModelName;
    }
 
    void writeAdjStates(
        const label writeMesh,
        const wordList& additionalOutput);
 
    scalar getElapsedClockTime()
    {
        return runTimePtr_->elapsedClockTime();
    }
 
    scalar getElapsedCpuTime()
    {
        return runTimePtr_->elapsedCpuTime();
    }
 
#ifdef CODI_ADR
 
    codi::RealReverse::Tape& globalADTape_;
 
#endif
 
    template<class classType>
    label validateField(const classType& field);
 
    template<class classType>
    label validateVectorField(const classType& field);
 
    label validateStates();
 
    void getInitStateVals(HashTable<scalar>& initState);
 
    void resetStateVals();
 
    void writeSensMapSurface(
        const word name,
        const double* dFdXs,
        const double* Xs,
        const label size,
        const double timeName);
 
    void writeSensMapField(
        const word name,
        const double* dFdField,
        const word fieldType,
        const double timeName);
 
    scalar getLatestTime()
    {
        instantList timeDirs = runTimePtr_->findTimes(runTimePtr_->path(), runTimePtr_->constant());
        scalar latestTime = timeDirs.last().value();
        return latestTime;
    }
 
    void writeAdjointFields(
        const word function,
        const double writeTime,
        const double* psi);
};
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
template<class classType>
label DASolver::validateField(const classType& field)
{
    /*
    Description:
        Check if a field variable has invalid in it. If the values are valid, return 0
    */
 
    forAll(field, idxI)
    {
        const scalar& val = field[idxI];
        if (std::isnan(val) || std::isinf(val) || fabs(val) > 1e15)
        {
            return 1;
        }
    }
    forAll(field.boundaryField(), patchI)
    {
        forAll(field.boundaryField()[patchI], faceI)
        {
            const scalar& val = field.boundaryField()[patchI][faceI];
            if (std::isnan(val) || std::isinf(val) || fabs(val) > 1e15)
            {
                return 1;
            }
        }
    }
 
    return 0;
}
 
template<class classType>
label DASolver::validateVectorField(const classType& field)
{
    /*
    Description:
        Check if a field variable has invalid in it. If the values are valid, return 0
    */
 
    forAll(field, idxI)
    {
        for (label compI = 0; compI < 3; compI++)
        {
            const scalar& val = field[idxI][compI];
            if (std::isnan(val) || std::isinf(val) || fabs(val) > 1e15)
            {
                return 1;
            }
        }
    }
    forAll(field.boundaryField(), patchI)
    {
        forAll(field.boundaryField()[patchI], faceI)
        {
            for (label compI = 0; compI < 3; compI++)
            {
                const scalar& val = field.boundaryField()[patchI][faceI][compI];
                if (std::isnan(val) || std::isinf(val) || fabs(val) > 1e15)
                {
                    return 1;
                }
            }
        }
    }
 
    return 0;
}
 
} // End namespace Foam
 
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
#endif
 
// ************************************************************************* //