1 /* 2 * This file is part of CasADi. 3 * 4 * CasADi -- A symbolic framework for dynamic optimization. 5 * Copyright (C) 2010-2014 Joel Andersson, Joris Gillis, Moritz Diehl, 6 * K.U. Leuven. All rights reserved. 7 * Copyright (C) 2011-2014 Greg Horn 8 * 9 * CasADi is free software; you can redistribute it and/or 10 * modify it under the terms of the GNU Lesser General Public 11 * License as published by the Free Software Foundation; either 12 * version 3 of the License, or (at your option) any later version. 13 * 14 * CasADi is distributed in the hope that it will be useful, 15 * but WITHOUT ANY WARRANTY; without even the implied warranty of 16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 17 * Lesser General Public License for more details. 18 * 19 * You should have received a copy of the GNU Lesser General Public 20 * License along with CasADi; if not, write to the Free Software 21 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA 22 * 23 */ 24 25 26 #ifndef CASADI_SUNDIALS_INTERFACE_HPP 27 #define CASADI_SUNDIALS_INTERFACE_HPP 28 29 #include <casadi/interfaces/sundials/casadi_sundials_common_export.h> 30 #include "casadi/core/integrator_impl.hpp" 31 32 #include <nvector/nvector_serial.h> 33 #include <sundials/sundials_dense.h> 34 #include <sundials/sundials_iterative.h> 35 #include <sundials/sundials_types.h> 36 37 #include <ctime> 38 39 /// \cond INTERNAL 40 namespace casadi { 41 42 // IdasMemory 43 struct CASADI_SUNDIALS_COMMON_EXPORT SundialsMemory : public IntegratorMemory { 44 // Current time 45 double t; 46 47 // N-vectors for the forward integration 48 N_Vector xz, xzdot, q; 49 50 // N-vectors for the backward integration 51 N_Vector rxz, rxzdot, rq; 52 53 // Initialize or reinitialize? 54 bool first_callB; 55 56 // Parameters 57 double *p, *rp; 58 59 // Jacobian 60 double *jac, *jacB; 61 62 /// Stats 63 long nsteps, nfevals, nlinsetups, netfails; 64 int qlast, qcur; 65 double hinused, hlast, hcur, tcur; 66 long nniters, nncfails; 67 68 long nstepsB, nfevalsB, nlinsetupsB, netfailsB; 69 int qlastB, qcurB; 70 double hinusedB, hlastB, hcurB, tcurB; 71 long nnitersB, nncfailsB; 72 73 // Temporaries for [x;z] or [rx;rz] 74 double *v1, *v2; 75 76 /// number of checkpoints stored so far 77 int ncheck; 78 79 /// Linear solver memory objects 80 int mem_linsolF, mem_linsolB; 81 82 /// Constructor 83 SundialsMemory(); 84 85 /// Destructor 86 ~SundialsMemory(); 87 }; 88 89 class CASADI_SUNDIALS_COMMON_EXPORT SundialsInterface : public Integrator { 90 public: 91 /** \brief Constructor */ 92 SundialsInterface(const std::string& name, const Function& dae); 93 94 /** \brief Destructor */ 95 ~SundialsInterface() override=0; 96 97 ///@{ 98 /** \brief Options */ 99 static const Options options_; get_options() const100 const Options& get_options() const override { return options_;} 101 ///@} 102 103 /** \brief Initialize */ 104 void init(const Dict& opts) override; 105 106 /** \brief Initalize memory block */ 107 int init_mem(void* mem) const override; 108 109 /** \brief Free memory block */ 110 void free_mem(void *mem) const override; 111 112 /** \brief Get relative tolerance */ get_reltol() const113 double get_reltol() const override { return reltol_;} 114 115 /** \brief Get absolute tolerance */ get_abstol() const116 double get_abstol() const override { return abstol_;} 117 118 // Get system Jacobian 119 virtual Function getJ(bool backward) const = 0; 120 121 /// Get all statistics 122 Dict get_stats(void* mem) const override; 123 124 /** \brief Print solver statistics */ 125 void print_stats(IntegratorMemory* mem) const override; 126 127 /** \brief Reset the forward problem and bring the time back to t0 */ 128 void reset(IntegratorMemory* mem, double t, const double* x, 129 const double* z, const double* p) const override; 130 131 /** \brief Reset the backward problem and take time to tf */ 132 void resetB(IntegratorMemory* mem, double t, const double* rx, 133 const double* rz, const double* rp) const override; 134 135 /** \brief Cast to memory object */ to_mem(void * mem)136 static SundialsMemory* to_mem(void *mem) { 137 SundialsMemory* m = static_cast<SundialsMemory*>(mem); 138 casadi_assert_dev(m); 139 return m; 140 } 141 142 ///@{ 143 /// Options 144 double abstol_, reltol_; 145 casadi_int max_num_steps_; 146 bool stop_at_end_; 147 bool quad_err_con_; 148 casadi_int steps_per_checkpoint_; 149 bool disable_internal_warnings_; 150 casadi_int max_multistep_order_; 151 std::string linear_solver_; 152 Dict linear_solver_options_; 153 casadi_int max_krylov_; 154 bool use_precon_; 155 bool second_order_correction_; 156 double step0_; 157 double max_step_size_; 158 double nonlin_conv_coeff_; 159 casadi_int max_order_; 160 ///@} 161 162 /// Linear solver 163 Linsol linsolF_, linsolB_; 164 165 /// Supported iterative solvers in Sundials 166 enum NewtonScheme {SD_DIRECT, SD_GMRES, SD_BCGSTAB, SD_TFQMR} newton_scheme_; 167 168 // Supported interpolations in Sundials 169 enum InterpType {SD_POLYNOMIAL, SD_HERMITE} interp_; 170 171 /// Linear solver data (dense) -- what is this? 172 struct LinSolDataDense {}; 173 174 /** \brief Set the (persistent) work vectors */ 175 void set_work(void* mem, const double**& arg, double**& res, 176 casadi_int*& iw, double*& w) const override; 177 178 // Print a variable printvar(const std::string & id,double v)179 static void printvar(const std::string& id, double v) { 180 uout() << id << " = " << v << std::endl; 181 } 182 // Print an N_Vector printvar(const std::string & id,N_Vector v)183 static void printvar(const std::string& id, N_Vector v) { 184 std::vector<double> tmp(NV_DATA_S(v), NV_DATA_S(v)+NV_LENGTH_S(v)); 185 uout() << id << " = " << tmp << std::endl; 186 } 187 188 /** \brief Serialize an object without type information */ 189 void serialize_body(SerializingStream &s) const override; 190 191 protected: 192 /** \brief Deserializing constructor */ 193 explicit SundialsInterface(DeserializingStream& s); 194 }; 195 196 // Check if N_Vector is regular is_regular(N_Vector v)197 inline bool is_regular(N_Vector v) { 198 std::vector<double> tmp(NV_DATA_S(v), NV_DATA_S(v)+NV_LENGTH_S(v)); 199 return is_regular(tmp); 200 } 201 202 } // namespace casadi 203 204 /// \endcond 205 #endif // CASADI_SUNDIALS_INTERFACE_HPP 206