Program listing for file kernel/src/simulationTools/Hem5OSI.hpp

Program listing for file kernel/src/simulationTools/Hem5OSI.hpp#

  1#ifndef Hem5OSI_H
  2#define Hem5OSI_H
  3
  4#include <vector>
  5
  6#include "OneStepIntegrator.hpp"
  7#include "SiconosConst.hpp"
  8#include "SiconosFortran.h"
  9
 10
 11class Hem5OSI : public OneStepIntegrator {
 12 private:
 13  ACCEPT_SERIALIZATION(Hem5OSI);
 14  static constexpr auto HEM5_ATOL_DEFAULT = 100 * siconos::internal::MACHINE_PREC;
 15  static constexpr auto HEM5_RTOL_DEFAULT = 10 * siconos::internal::MACHINE_PREC;
 16  static constexpr double INITIAL_GUESS_TS = 1.e-3;
 17
 18
 19  std::vector<int> _intData = {0, 0, 0, 0, 0, 0, 0, 0, 0};
 20
 21  int _idid{0};
 22
 23
 24  std::vector<double> rtol = {};
 25
 26  std::vector<double> atol = {};
 27
 28  std::vector<double> rwork = {};
 29
 30  std::vector<int> iwork = {};
 31
 32  double _timeStep{INITIAL_GUESS_TS};
 33
 34
 35  SP::BlockVector _qWork;
 36
 37  SP::BlockVector _vWork;
 38
 39  SP::BlockVector _uWork;
 40
 41  SP::BlockVector _aWork;
 42
 43  SP::BlockVector _lambdaWork;
 44
 45  SP::BlockVector _forcesWork;
 46
 47  SP::SiconosVector _qtmp;
 48  SP::SiconosVector _vtmp;
 49  SP::SiconosVector _utmp;
 50  SP::SiconosVector _atmp;
 51  SP::SiconosVector _lambdatmp;
 52  SP::SiconosVector _forcestmp;
 53
 54
 55  struct _NSLEffectOnFreeOutput;
 56  friend struct _NSLEffectOnFreeOutput;
 57
 58  class Hem5OSI_impl {
 59   public:
 60    std::shared_ptr<Hem5OSI> hem5osi{nullptr};
 61    void fprob(int* IFCN, int* NQ, int* NV, int* NU, int* NL, int* LDG, int* LDF, int* LDA,
 62               int* NBLK, int* NMRC, int* NPGP, int* NPFL, int* INDGR, int* INDGC, int* INDFLR,
 63               int* INDFLC, double* time, double* q, double* v, double* u, double* xl,
 64               double* G, double* GQ, double* F, double* GQQ, double* GT, double* FL,
 65               double* QDOT, double* UDOT, double* AM);
 66    void solout(int* MODE, int* NSTEP, int* NQ, int* NV, int* NU, int* NL, int* LDG, int* LDF,
 67                int* LDA, int* LRDO, int* LIDO, siconos::fortran::hairer::fprobpointer FPROB,
 68                double* q, double* v, double* u, double* DOWK, int* IDOWK);
 69
 70    Hem5OSI_impl(std::shared_ptr<Hem5OSI> h) : hem5osi(h) {}
 71  };
 72
 73 public:
 74  std::shared_ptr<Hem5OSI_impl> _impl{nullptr};
 75
 76  enum Hem5OSI_ds_workVector_id { FREE, WORK_LENGTH };
 77
 78  enum Hem5OSI_interaction_workVector_id { OSNSP_RHS, WORK_INTERACTION_LENGTH };
 79
 80  enum Hem5OSI_interaction_workBlockVector_id { xfree, BLOCK_WORK_LENGTH };
 81
 82
 83  Hem5OSI();
 84
 85
 86  ~Hem5OSI() noexcept = default;
 87
 88
 89  inline const std::vector<int> intData() const { return _intData; }
 90
 91
 92  inline int intData(unsigned int i) const { return _intData[i]; }
 93
 94  inline void setIntData(unsigned int i, int newValue) { _intData[i] = newValue; }
 95
 96
 97  inline const std::vector<double>& getRtol() const { return rtol; }
 98
 99
100  inline const std::vector<double>& getAtol() const { return atol; }
101
102
103  inline int getMaxNstep() const { return iwork[11]; }
104
105
106  inline const std::vector<double>& getRwork() const { return rwork; }
107
108
109  inline const std::vector<int>& getIwork() const { return iwork; }
110
111
112  void setTol(int itol, std::vector<double>&& rtol, std::vector<double>&& atol);
113
114
115  void setTol(int itol, double rtol, double atol);
116
117
118  void setMaxNstep(int nstepmax);
119
120
121  void setMaxStepSize(double maxstepsize);
122
123
124  void updateIntData();
125
126
127  void updateData();
128
129
130  void fillqWork(int* sizex, double* x);
131
132
133  void fillvWork(int* sizex, double* x);
134
135
136  void computeRhs(double);
137
138
139  void computeJacobianRhs(double);
140
141  unsigned int numberOfConstraints();
142
143  void f(int* sizeOfX, double* time, double* x, double* xdot);
144
145  void g(int* nEq, double* time, double* x, int* ng, double* gOut);
146
147  void jacobianfx(int*, double*, double*, int*, int*, double*, int*);
148
149
150  void initialize() override;
151
152  void initializeWorkVectorsForDS(double t, SP::DynamicalSystem ds) override;
153
154
155  void initializeWorkVectorsForInteraction(Interaction& inter,
156                                           InteractionProperties& interProp,
157                                           DynamicalSystemsGraph& DSG) override;
158
159
160  unsigned int numberOfIndexSets() const override { return 3; };
161
162
163  void integrate(double& tinit, double& tend, double& tout, int& idid) override;
164
165
166  void updateState(const unsigned int level) override;
167
168  void prepareNewtonIteration(double time) override { assert(0); };
169
170
171  void computeFreeOutput(InteractionsGraph::VDescriptor& vertex_inter,
172                         OneStepNSProblem* osnsp) override;
173
174
175  SiconosVector& osnsp_rhs(InteractionsGraph::VDescriptor& vertex_inter,
176                           InteractionsGraph& indexSet) override {
177    return *(*indexSet.properties(vertex_inter).workVectors)[Hem5OSI::OSNSP_RHS];
178  };
179
180
181  void display() override;
182
183  ACCEPT_STD_VISITORS();
184};
185
186#endif