ODE.h

Go to the documentation of this file.

#ifndef ODE_H_
#define ODE_H_
 
#include <boost/property_tree/ptree.hpp>
 
#include <cvodes/cvodes.h>               /* access to CVODE                 */
#include <nvector/nvector_serial.h>    /* access to serial N_Vector       */
#include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix       */
#include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver */
#include <sunnonlinsol/sunnonlinsol_newton.h>     /* access to the newton SUNNonlinearSolver      */
#include <sunnonlinsol/sunnonlinsol_fixedpoint.h> /* access to the fixed point SUNNonlinearSolver */
 
 
#include "MUQ/Modeling/ModPiece.h"
 
namespace muq {
namespace Modeling {
 
 
  class ODE : public ModPiece {
  public:
 
 
    ODE(std::shared_ptr<ModPiece>   const& rhsIn,
        Eigen::VectorXd             const& evalTimesIn,
        boost::property_tree::ptree        opts = boost::property_tree::ptree()) : ODE(rhsIn, evalTimesIn, true, opts){};
 
    ODE(std::shared_ptr<ModPiece>   const& rhsIn,
        Eigen::VectorXd             const& evalTimesIn,
        bool                               isAuto,
        boost::property_tree::ptree        opts = boost::property_tree::ptree());
 
 
    virtual ~ODE() = default;
 
    const bool isAutonomous;
 
    const double startTime;
 
    protected:
 
      virtual void EvaluateImpl(ref_vector<Eigen::VectorXd> const& inputs) override;
 
      virtual void GradientImpl(unsigned int                       outWrt,
                                unsigned int                       inWrt,
                                ref_vector<Eigen::VectorXd> const& inputs,
                                Eigen::VectorXd const&             sens) override;
 
      virtual void JacobianImpl(unsigned int outWrt,
                                unsigned int inWrt,
                                ref_vector<Eigen::VectorXd> const& inputs) override;
 
      virtual void ApplyJacobianImpl(unsigned int outWrt,
                                      unsigned int inWrt,
                                      ref_vector<Eigen::VectorXd> const& inputs,
                                      Eigen::VectorXd const& vec) override;
 
      struct NonlinearSolverOptions{
        NonlinearSolverOptions() = default;
        NonlinearSolverOptions(boost::property_tree::ptree const& opts);
 
        void SetOptions(boost::property_tree::ptree const& opts);
 
        std::string method; // e.g., Iter, Newton
        int maxIters;
      };
 
      struct LinearSolverOptions{
        LinearSolverOptions() = default;
        LinearSolverOptions(boost::property_tree::ptree const& opts);
        void SetOptions(boost::property_tree::ptree const& opts);
 
        std::string method; // e.g., Dense, SPGMR, SPFGMR, SPBCGS, SPTFQMR, PCG
        int maxIters; // maximum number of iterations allowed for iterative solvers
 
      };
 
      struct IntegratorOptions{
        IntegratorOptions() = default;
        IntegratorOptions(boost::property_tree::ptree const& opts);
 
        void SetOptions(boost::property_tree::ptree const& opts);
 
        std::string method; // e.g., Adams or BDF
        double reltol; // relative tolerance
        double abstol; // absolute tolerance
        double maxStepSize; // maximum allowed step inputSize
        int maxNumSteps; // maximum allowed steps
        int checkPtGap; // gap between check points
 
        double sensRelTol; // relative tolerance for forward sensitivities
        double sensAbsTol; // absolute tolerance for forward sensitivities.
 
        double adjRelTol; // relative tolerance for adjoint sensitivities.
        double adjAbsTol; // absolute tolerance for adjoint sensitivities.
      };
 
      struct Interface{
 
        static int RHS(realtype time, N_Vector state, N_Vector deriv, void *user_data);
        static int RHSJacobian(realtype time, N_Vector state, N_Vector rhs, SUNMatrix jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
        static int RHSJacobianAction(N_Vector v, N_Vector Jv, realtype time, N_Vector state, N_Vector rhs, void *user_data, N_Vector tmp);
 
        static int SensitivityRHS(int Ns, realtype time, N_Vector y, N_Vector ydot, N_Vector *ys, N_Vector *ySdot, void *user_data, N_Vector tmp1, N_Vector tmp2);
 
        static int AdjointRHS(double time, N_Vector state, N_Vector lambda, N_Vector deriv, void *user_data);
        static int AdjointJacobian(double time, N_Vector state, N_Vector lambda, N_Vector rhs, SUNMatrix jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
 
        static int AdjointJacobianAction(N_Vector target, N_Vector output, double time, N_Vector state, N_Vector lambda, N_Vector adjRhs, void *user_data, N_Vector tmp);
 
        static int AdjointQuad(realtype time, N_Vector state, N_Vector lambda, N_Vector quadRhs, void *user_data);
      };
 
      NonlinearSolverOptions nonlinOpts;
      LinearSolverOptions linOpts;
      IntegratorOptions intOpts;
 
      // Store the times when where we want to return the solution of the ODE
      const Eigen::VectorXd evalTimes;
 
      // Store a ModPiece that evaluates the RHS function
      std::shared_ptr<ModPiece>  rhs;
 
    private:
 
      std::pair<SUNLinearSolver, SUNMatrix> CreateLinearSolver(void *cvode_mem, N_Vector const& state);
 
      SUNNonlinearSolver CreateNonlinearSolver(void *cvode_mem, N_Vector const& state);
 
      std::pair<SUNLinearSolver, SUNMatrix> CreateAdjointLinearSolver(void *cvode_mem, N_Vector const& state, int indexB);
 
      SUNNonlinearSolver CreateAdjointNonlinearSolver(void *cvode_mem, N_Vector const& state, int indexB);
 
 
      static Eigen::VectorXi GetInputSizes(std::shared_ptr<ModPiece> const& rhsIn,
                                           bool                             isAuto);
 
      static Eigen::VectorXi GetOutputSizes(std::shared_ptr<ModPiece> const& rhsIn,
                                            Eigen::VectorXd           const& evalTimes);
 
      static double GetStartTime(Eigen::VectorXd             const& evalTimesIn,
                                 boost::property_tree::ptree const& opts);
 
      static void CheckFlag(void *returnvalue, const char *funcname, int opt);
    };
 
} // namespace Modeling
} // namespace muq
 
#endif // #ifndef ODE_H_