idas_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_IDAS_INTERFACE_HPP
#define CASADI_IDAS_INTERFACE_HPP
 
#include <casadi/interfaces/sundials/casadi_integrator_idas_export.h>
#include "sundials_interface.hpp"
#include <idas/idas.h>            /* prototypes for CVODE fcts. and consts. */
#include <idas/idas_dense.h>
#include <idas/idas_band.h>
#include <idas/idas_spgmr.h>
#include <idas/idas_spbcgs.h>
#include <idas/idas_sptfqmr.h>
#include <idas/idas_impl.h> /* Needed for the provided linear solver */
#include <ctime>
 
 
namespace casadi {
 
// Forward declaration
class IdasInterface;
 
// IdasMemory
struct  IdasMemory : public SundialsMemory {
  const IdasInterface& self;
 
  void* mem;
 
  int whichB;
 
  double cj_last;
 
  IdasMemory(const IdasInterface& s);
 
  ~IdasMemory();
};
 
class  IdasInterface : public SundialsInterface {
 public:
 
  IdasInterface(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 IdasInterface(name, dae, t0, tout);
  }
 
  ~IdasInterface() override;
 
  // Get name of the plugin
  const char* plugin_name() const override { return "idas";}
 
  // Get name of the class
  std::string class_name() const override { return "IdasInterface";}
 
 
  static const Options options_;
  const Options& get_options() const override { return options_;}
 
  void init(const Dict& opts) override;
 
  void* alloc_mem() const override { return new IdasMemory(*this);}
 
  int init_mem(void* mem) const override;
 
  void free_mem(void *mem) const override { delete static_cast<IdasMemory*>(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 resetB(IntegratorMemory* mem) const override;
 
  void z_impulseB(IdasMemory* m, const double* rz) const;
 
  int solve_transposed(IdasMemory* m, double t, const double* xz, const double* rxz,
    const double* rhs, double* sol) const;
 
  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 IdasMemory* to_mem(void *mem) {
    IdasMemory* m = static_cast<IdasMemory*>(mem);
    casadi_assert_dev(m);
    return m;
  }
 
  static const std::string meta_doc;
 
 protected:
 
  // Sundials callback functions
  static int resF(double t, N_Vector xz, N_Vector xzdot, N_Vector rr, void *user_data);
  static int resB(double t, N_Vector xz, N_Vector xzdot, N_Vector rxz, N_Vector rxzdot,
    N_Vector rr, void *user_data);
  static void ehfun(int error_code, const char *module, const char *function, char *msg,
    void *eh_data);
  static int jtimesF(double t, N_Vector xz, N_Vector xzdot, N_Vector rr, N_Vector v,
    N_Vector Jv, double cj, void *user_data, N_Vector tmp1, N_Vector tmp2);
  static int jtimesB(double t, N_Vector xz, N_Vector xzdot, N_Vector xzB, N_Vector xzdotB,
    N_Vector resvalB, N_Vector vB, N_Vector JvB, double cjB,
    void *user_data, N_Vector tmp1B, N_Vector tmp2B);
  static int rhsQF(double t, N_Vector xz, N_Vector xzdot, N_Vector qdot, void *user_data);
  static int rhsQB(double t, N_Vector xz, N_Vector xzdot, N_Vector rxz,
    N_Vector rxzdot, N_Vector ruqdot, void *user_data);
  static int psolveF(double t, N_Vector xz, N_Vector xzdot, N_Vector rr, N_Vector rvec,
    N_Vector zvec, double cj, double delta, void *user_data, N_Vector tmp);
  static int psetupF(double t, N_Vector xz, N_Vector xzdot, N_Vector rr, double cj,
    void* user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
  static int psolveB(double t, N_Vector xz, N_Vector xzdot, N_Vector rxz, N_Vector rxzdot,
    N_Vector resvalB, N_Vector rvecB, N_Vector zvecB, double cjB,
    double deltaB, void *user_dataB, N_Vector tmpB);
  static int psetupB(double t, N_Vector xz, N_Vector xzdot, N_Vector rxz, N_Vector rxzdot,
    N_Vector resvalB, double cjB, void *user_dataB,
    N_Vector tmp1B, N_Vector tmp2B, N_Vector tmp3B);
  static int lsetupF(IDAMem IDA_mem, N_Vector xz, N_Vector xzdot, N_Vector resp,
    N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3);
  static int lsolveF(IDAMem IDA_mem, N_Vector b, N_Vector weight, N_Vector ycur,
    N_Vector xzdotcur, N_Vector rescur);
  static int lsetupB(IDAMem IDA_mem, N_Vector xz, N_Vector xzdot, N_Vector resp,
    N_Vector vtemp1, N_Vector vtemp2, N_Vector vtemp3);
  static int lsolveB(IDAMem IDA_mem, N_Vector b, N_Vector weight, N_Vector ycur,
    N_Vector xzdotcur, N_Vector rescur);
 public:
 
  // Throw error
  static void idas_error(const char* module, int flag);
 
  // Options
  bool cj_scaling_;
  bool calc_ic_;
  bool calc_icB_;
  bool suppress_algebraic_;
  std::vector<double> abstolv_;
  double first_time_;
 
  // Constraints on decision variables
  std::vector<casadi_int> y_c_;
 
  //  Initial values for \p xdot
  std::vector<double> init_xdot_;
 
  void serialize_body(SerializingStream &s) const override;
 
  static ProtoFunction* deserialize(DeserializingStream& s) { return new IdasInterface(s); }
 
 protected:
  explicit IdasInterface(DeserializingStream& s);
};
 
} // namespace casadi
 
#endif // CASADI_IDAS_INTERFACE_HPP