cvodes_interface.hpp

/*
 *    This file is part of CasADi.
 *
 *    CasADi -- A symbolic framework for dynamic optimization.
 *    Copyright (C) 2010-2023 Joel Andersson, Joris Gillis, Moritz Diehl,
 *                            KU Leuven. All rights reserved.
 *    Copyright (C) 2011-2014 Greg Horn
 *
 *    CasADi is free software; you can redistribute it and/or
 *    modify it under the terms of the GNU Lesser General Public
 *    License as published by the Free Software Foundation; either
 *    version 3 of the License, or (at your option) any later version.
 *
 *    CasADi is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 *    Lesser General Public License for more details.
 *
 *    You should have received a copy of the GNU Lesser General Public
 *    License along with CasADi; if not, write to the Free Software
 *    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
 *
 */
 
 
#ifndef CASADI_CVODES_INTERFACE_HPP
#define CASADI_CVODES_INTERFACE_HPP
 
#include <casadi/interfaces/sundials/casadi_integrator_cvodes_export.h>
#include "sundials_interface.hpp"
#include <cvodes/cvodes.h>            /* prototypes for CVode fcts. and consts. */
#include <cvodes/cvodes_dense.h>
#include <cvodes/cvodes_band.h>
#include <cvodes/cvodes_spgmr.h>
#include <cvodes/cvodes_spbcgs.h>
#include <cvodes/cvodes_sptfqmr.h>
#include <cvodes/cvodes_impl.h> /* Needed for the provided linear solver */
#include <ctime>
 
namespace casadi {
 
// Forward declaration
class CvodesInterface;
 
// CvodesMemory
struct  CvodesMemory : public SundialsMemory {
  const CvodesInterface& self;
 
  // CVodes memory block
  void* mem;
 
  // Ids of backward problem
  int whichB;
 
  // Remember the gamma and gammaB from last factorization
  double gamma, gammaB;
 
  CvodesMemory(const CvodesInterface& s);
 
  ~CvodesMemory();
};
 
class  CvodesInterface : public SundialsInterface {
 public:
  CvodesInterface(const std::string& name, const Function& dae, double t0,
    const std::vector<double>& tout);
 
  static Integrator* creator(const std::string& name, const Function& dae,
      double t0, const std::vector<double>& tout) {
    return new CvodesInterface(name, dae, t0, tout);
  }
 
  ~CvodesInterface() override;
 
  // Get name of the plugin
  const char* plugin_name() const override { return "cvodes";}
 
  // Get name of the class
  std::string class_name() const override { return "CvodesInterface";}
 
 
  static const Options options_;
  const Options& get_options() const override { return options_;}
  double min_step_size_;
  bool always_recalculate_jacobian_;
 
 
  void init(const Dict& opts) override;
 
  void* alloc_mem() const override { return new CvodesMemory(*this);}
 
  int init_mem(void* mem) const override;
 
  void free_mem(void *mem) const override { delete static_cast<CvodesMemory*>(mem);}
 
  void reset(IntegratorMemory* mem,
    const double* u, const double* x, const double* z, const double* p) const override;
 
  void advance(IntegratorMemory* mem,
    const double* u, double* x, double* z, double* q) const override;
 
  void impulseB(IntegratorMemory* mem,
    const double* rx, const double* rz, const double* rp) const override;
 
  void retreat(IntegratorMemory* mem, const double* u,
    double* rx, double* rq, double* uq) const override;
 
  static CvodesMemory* to_mem(void *mem) {
    CvodesMemory* m = static_cast<CvodesMemory*>(mem);
    casadi_assert_dev(m);
    return m;
  }
 
  static const std::string meta_doc;
 
 protected:
 
  // Sundials callback functions
  static void ehfun(int error_code, const char *module, const char *function, char *msg,
    void *user_data);
  static int rhsF(double t, N_Vector x, N_Vector xdot, void *user_data);
  static int rhsB(double t, N_Vector x, N_Vector xB, N_Vector xdotB, void *user_data);
  static int rhsQF(double t, N_Vector x, N_Vector qdot, void *user_data);
  static int rhsQB(double t, N_Vector x, N_Vector rx, N_Vector ruqdot, void *user_data);
  static int jtimesF(N_Vector v, N_Vector Jv, double t, N_Vector x, N_Vector xdot,
    void *user_data, N_Vector tmp);
  static int jtimesB(N_Vector vB, N_Vector JvB, double t, N_Vector x, N_Vector xB,
    N_Vector xdotB, void *user_data , N_Vector tmpB);
  static int psolveF(double t, N_Vector x, N_Vector xdot, N_Vector r, N_Vector z,
    double gamma, double delta, int lr, void *user_data, N_Vector tmp);
  static int psolveB(double t, N_Vector x, N_Vector xB, N_Vector xdotB, N_Vector rvecB,
    N_Vector zvecB, double gammaB, double deltaB,
    int lr, void *user_data, N_Vector tmpB);
  static int psetupF(double t, N_Vector x, N_Vector xdot, booleantype jok,
    booleantype *jcurPtr, double gamma, void *user_data,
    N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
  static int psetupB(double t, N_Vector x, N_Vector xB, N_Vector xdotB,
    booleantype jokB, booleantype *jcurPtrB, double gammaB,
    void *user_data, N_Vector tmp1B, N_Vector tmp2B, N_Vector tmp3B);
  static int lsetupF(CVodeMem cv_mem, int convfail, N_Vector x, N_Vector xdot,
    booleantype *jcurPtr, N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3);
  static int lsolveF(CVodeMem cv_mem, N_Vector b, N_Vector weight, N_Vector x, N_Vector xdot);
  static int lsetupB(CVodeMem cv_mem, int convfail, N_Vector x, N_Vector xdot,
    booleantype *jcurPtr, N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3);
  static int lsolveB(CVodeMem cv_mem, N_Vector b, N_Vector weight, N_Vector x, N_Vector xdot);
 
  // Throw error
  static void cvodes_error(const char* module, int flag);
 
  casadi_int lmm_; // linear multistep method
  casadi_int iter_; // nonlinear solver iteration
 
 public:
 
  void serialize_body(SerializingStream &s) const override;
 
  static ProtoFunction* deserialize(DeserializingStream& s) { return new CvodesInterface(s); }
 
 protected:
 
  explicit CvodesInterface(DeserializingStream& s);
};
 
} // namespace casadi
 
#endif // CASADI_CVODES_INTERFACE_HPP