Simbody  3.5
SimTKcpodes.h
Go to the documentation of this file.
1 #ifndef SimTK_SIMMATH_CPODES_H_
2 #define SimTK_SIMMATH_CPODES_H_
3 
4 /* -------------------------------------------------------------------------- *
5  * Simbody(tm): SimTKmath *
6  * -------------------------------------------------------------------------- *
7  * This is part of the SimTK biosimulation toolkit originating from *
8  * Simbios, the NIH National Center for Physics-Based Simulation of *
9  * Biological Structures at Stanford, funded under the NIH Roadmap for *
10  * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. *
11  * *
12  * Portions copyright (c) 2006-12 Stanford University and the Authors. *
13  * Authors: Michael Sherman *
14  * Contributors: *
15  * *
16  * Licensed under the Apache License, Version 2.0 (the "License"); you may *
17  * not use this file except in compliance with the License. You may obtain a *
18  * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
19  * *
20  * Unless required by applicable law or agreed to in writing, software *
21  * distributed under the License is distributed on an "AS IS" BASIS, *
22  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
23  * See the License for the specific language governing permissions and *
24  * limitations under the License. *
25  * -------------------------------------------------------------------------- */
26 
33 #include "SimTKcommon.h"
35 
36 #include <cstdio> // Needed for "FILE".
37 
38 namespace SimTK {
39 
50 public:
51  virtual ~CPodesSystem() {}
52  // The default implementations of these virtual functions
53  // just throw an "unimplemented virtual function" exception.
54 
55  // At least one of these two must be supplied by the concrete class.
56  virtual int explicitODE(Real t, const Vector& y, Vector& fout) const;
57  virtual int implicitODE(Real t, const Vector& y, const Vector& yp,
58  Vector& fout) const;
59 
60  virtual int constraint(Real t, const Vector& y,
61  Vector& cout) const;
62  virtual int project(Real t, const Vector& ycur, Vector& corr,
63  Real epsProj, Vector& err) const; // err is in/out
64  virtual int quadrature(Real t, const Vector& y,
65  Vector& qout) const;
66  virtual int root(Real t, const Vector& y, const Vector& yp,
67  Vector& gout) const;
68  virtual int weight(const Vector& y, Vector& weights) const;
69  virtual void errorHandler(int error_code, const char* module,
70  const char* function, char* msg) const;
71 
72  //TODO: Jacobian functions
73 };
74 
75 
76 // These static functions are private to the current (client-side) compilation
77 // unit. They are used to navigate the client-side CPodesSystem virtual function
78 // table, which cannot be done on the library side. Note that these are defined
79 // in the SimTK namespace so don't need "SimTK" in their names.
80 static int explicitODE_static(const CPodesSystem& sys,
81  Real t, const Vector& y, Vector& fout)
82  { return sys.explicitODE(t,y,fout); }
83 
84 static int implicitODE_static(const CPodesSystem& sys,
85  Real t, const Vector& y, const Vector& yp, Vector& fout)
86  { return sys.implicitODE(t,y,yp,fout); }
87 
88 static int constraint_static(const CPodesSystem& sys,
89  Real t, const Vector& y, Vector& cout)
90  { return sys.constraint(t,y,cout); }
91 
92 static int project_static(const CPodesSystem& sys,
93  Real t, const Vector& ycur, Vector& corr,
94  Real epsProj, Vector& err)
95  { return sys.project(t,ycur,corr,epsProj,err); }
96 
97 static int quadrature_static(const CPodesSystem& sys,
98  Real t, const Vector& y, Vector& qout)
99  { return sys.quadrature(t,y,qout); }
100 
101 static int root_static(const CPodesSystem& sys,
102  Real t, const Vector& y, const Vector& yp, Vector& gout)
103  { return sys.root(t,y,yp,gout); }
104 
105 static int weight_static(const CPodesSystem& sys,
106  const Vector& y, Vector& weights)
107  { return sys.weight(y,weights); }
108 
109 static void errorHandler_static(const CPodesSystem& sys,
110  int error_code, const char* module,
111  const char* function, char* msg)
112  { sys.errorHandler(error_code,module,function,msg); }
113 
123 public:
124  // no default constructor
125  // copy constructor and default assignment are suppressed
126 
127  enum ODEType {
128  UnspecifiedODEType=0,
130  ImplicitODE
131  };
132 
134  UnspecifiedLinearMultistepMethod=0,
136  Adams
137  };
138 
140  UnspecifiedNonlinearSystemIterationType=0,
142  Functional
143  };
144 
146  UnspecifiedToleranceType=0,
149  WeightFunction
150  };
151 
153  UnspecifiedProjectionNorm=0,
155  ErrorNorm
156  };
157 
159  UnspecifiedConstraintLinearity=0,
161  Nonlinear
162  };
163 
165  UnspecifiedProjectionFactorizationType=0,
169  ProjectWithQRPivot // for handling redundancy
170  };
171 
172  enum StepMode {
173  UnspecifiedStepMode=0,
177  OneStepTstop
178  };
179 
180  explicit CPodes
181  (ODEType ode=UnspecifiedODEType,
182  LinearMultistepMethod lmm=UnspecifiedLinearMultistepMethod,
183  NonlinearSystemIterationType nls=UnspecifiedNonlinearSystemIterationType)
184  {
185  // Perform construction of the CPodesRep on the library side.
186  librarySideCPodesConstructor(ode, lmm, nls);
187  // But fill in function pointers from the client side.
188  clientSideCPodesConstructor();
189  }
190 
191  // Values for 'flag' return values. These are just the "normal" return
192  // values; there are many more which are all negative and represent
193  // error conditions.
194  static const int Success = 0;
195  static const int TstopReturn = 1;
196  static const int RootReturn = 2;
197  static const int Warning = 99;
198  static const int TooMuchWork = -1;
199  static const int TooClose = -27;
200 
201  // These values should be used by user routines. "Success" is the
202  // same as above. A positive return value means "recoverable error",
203  // i.e., CPodes should cut the step size and try again, while a
204  // negative number means "unrecoverable error" which will kill
205  // CPODES altogether with a CP_xxx_FAIL error. The particular numerical
206  // values here have no significance, just + vs. -.
207  static const int RecoverableError = 9999;
208  static const int UnrecoverableError = -9999;
209 
210  ~CPodes();
211 
212  // Depending on the setting of ode_type at construction, init()
213  // and reInit() will tell CPodes to use either the explicitODE()
214  // or implicitODE() function from the CPodesSystem, so the user
215  // MUST have overridden at least one of those virtual methods.
216  int init(CPodesSystem& sys,
217  Real t0, const Vector& y0, const Vector& yp0,
218  ToleranceType tt, Real reltol, void* abstol);
219 
220  int reInit(CPodesSystem& sys,
221  Real t0, const Vector& y0, const Vector& yp0,
222  ToleranceType tt, Real reltol, void* abstol);
223 
224  // This tells CPodes to make use of the user's constraint()
225  // method from CPodesSystem, and perform projection internally.
226  int projInit(ProjectionNorm, ConstraintLinearity,
227  const Vector& ctol);
228 
229  // This tells CPodes to make use of the user's project()
230  // method from CPodesSystem.
231  int projDefine();
232 
233  // These tell CPodes to make use of the user's quadrature()
234  // method from CPodesSystem.
235  int quadInit(const Vector& q0);
236  int quadReInit(const Vector& q0);
237 
238  // This tells CPodes to make use of the user's root() method
239  // from CPodesSystem.
240  int rootInit(int nrtfn);
241 
242  // This tells CPodes to make use of the user's errorHandler()
243  // method from CPodesSystem.
244  int setErrHandlerFn();
245 
246  // These tells CPodes to make use of the user's weight()
247  // method from CPodesSystem.
248  int setEwtFn();
249 
250  // TODO: these routines should enable methods that are defined
251  // in the CPodesSystem, but a proper interface to the Jacobian
252  // routines hasn't been implemented yet.
253  int dlsSetJacFn(void* jac, void* jac_data);
254  int dlsProjSetJacFn(void* jacP, void* jacP_data);
255 
256 
257  int step(Real tout, Real* tret,
258  Vector& y_inout, Vector& yp_inout, StepMode=Normal);
259 
260  int setErrFile(FILE* errfp);
261  int setMaxOrd(int maxord);
262  int setMaxNumSteps(int mxsteps);
263  int setMaxHnilWarns(int mxhnil);
264  int setStabLimDet(bool stldet) ;
265  int setInitStep(Real hin);
266  int setMinStep(Real hmin);
267  int setMaxStep(Real hmax);
268  int setStopTime(Real tstop);
269  int setMaxErrTestFails(int maxnef);
270 
271  int setMaxNonlinIters(int maxcor);
272  int setMaxConvFails(int maxncf);
273  int setNonlinConvCoef(Real nlscoef);
274 
275  int setProjUpdateErrEst(bool proj_err);
276  int setProjFrequency(int proj_freq);
277  int setProjTestCnstr(bool test_cnstr);
278  int setProjLsetupFreq(int proj_lset_freq);
279  int setProjNonlinConvCoef(Real prjcoef);
280 
281  int setQuadErrCon(bool errconQ,
282  int tol_typeQ, Real reltolQ, void* abstolQ);
283 
284  int setTolerances(int tol_type, Real reltol, void* abstol);
285 
286  int setRootDirection(Array_<int>& rootdir);
287 
288  int getDky(Real t, int k, Vector& dky);
289 
290  int getQuad(Real t, Vector& yQout);
291  int getQuadDky(Real t, int k, Vector& dky);
292 
293  int getWorkSpace(int* lenrw, int* leniw);
294  int getNumSteps(int* nsteps);
295  int getNumFctEvals(int* nfevals);
296  int getNumLinSolvSetups(int* nlinsetups);
297  int getNumErrTestFails(int* netfails);
298  int getLastOrder(int* qlast);
299  int getCurrentOrder(int* qcur);
300  int getNumStabLimOrderReds(int* nslred);
301  int getActualInitStep(Real* hinused);
302  int getLastStep(Real* hlast);
303  int getCurrentStep(Real* hcur);
304  int getCurrentTime(Real* tcur);
305  int getTolScaleFactor(Real* tolsfac);
306  int getErrWeights(Vector& eweight);
307  int getEstLocalErrors(Vector& ele) ;
308  int getNumGEvals(int* ngevals);
309  int getRootInfo(int* rootsfound);
310  int getRootWindow(Real* tLo, Real* tHi);
311  int getIntegratorStats(int* nsteps,
312  int* nfevals, int* nlinsetups,
313  int* netfails, int* qlast,
314  int* qcur, Real* hinused, Real* hlast,
315  Real* hcur, Real* tcur);
316 
317  int getNumNonlinSolvIters(int* nniters);
318  int getNumNonlinSolvConvFails(int* nncfails);
319  int getNonlinSolvStats(int* nniters, int* nncfails);
320  int getProjNumProj(int* nproj);
321  int getProjNumCnstrEvals(int* nce);
322  int getProjNumLinSolvSetups(int* nsetupsP);
323  int getProjNumFailures(int* nprf) ;
324  int getProjStats(int* nproj,
325  int* nce, int* nsetupsP,
326  int* nprf);
327  int getQuadNumFunEvals(int* nqevals);
328  int getQuadErrWeights(Vector& eQweight);
329  char* getReturnFlagName(int flag);
330 
331 
332  int dlsGetWorkSpace(int* lenrwLS, int* leniwLS);
333  int dlsGetNumJacEvals(int* njevals);
334  int dlsGetNumFctEvals(int* nfevalsLS);
335  int dlsGetLastFlag(int* flag);
336  char* dlsGetReturnFlagName(int flag);
337 
338  int dlsProjGetNumJacEvals(int* njPevals);
339  int dlsProjGetNumFctEvals(int* ncevalsLS);
340 
341  int lapackDense(int N);
342  int lapackBand(int N, int mupper, int mlower);
343  int lapackDenseProj(int Nc, int Ny, ProjectionFactorizationType);
344 
345 private:
346  // This is how we get the client-side virtual functions to
347  // be callable from library-side code while maintaining binary
348  // compatibility.
349  typedef int (*ExplicitODEFunc)(const CPodesSystem&,
350  Real t, const Vector& y, Vector& fout);
351  typedef int (*ImplicitODEFunc)(const CPodesSystem&,
352  Real t, const Vector& y, const Vector& yp,
353  Vector& fout);
354  typedef int (*ConstraintFunc) (const CPodesSystem&,
355  Real t, const Vector& y, Vector& cout);
356  typedef int (*ProjectFunc) (const CPodesSystem&,
357  Real t, const Vector& ycur, Vector& corr,
358  Real epsProj, Vector& err);
359  typedef int (*QuadratureFunc) (const CPodesSystem&,
360  Real t, const Vector& y, Vector& qout);
361  typedef int (*RootFunc) (const CPodesSystem&,
362  Real t, const Vector& y, const Vector& yp,
363  Vector& gout);
364  typedef int (*WeightFunc) (const CPodesSystem&,
365  const Vector& y, Vector& weights);
366  typedef void (*ErrorHandlerFunc)(const CPodesSystem&,
367  int error_code, const char* module,
368  const char* function, char* msg);
369 
370  // Note that these routines do not tell CPodes to use the supplied
371  // functions. They merely provide the client-side addresses of functions
372  // which understand how to find the user's virtual functions, should those
373  // actually be provided. Control over whether to actually call any of these
374  // is handled elsewhere, with user-visible methods. These private methods
375  // are to be called only upon construction of the CPodes object here. They
376  // are not even dependent on which user-supplied concrete CPodesSystem is
377  // being used.
378  void registerExplicitODEFunc(ExplicitODEFunc);
379  void registerImplicitODEFunc(ImplicitODEFunc);
380  void registerConstraintFunc(ConstraintFunc);
381  void registerProjectFunc(ProjectFunc);
382  void registerQuadratureFunc(QuadratureFunc);
383  void registerRootFunc(RootFunc);
384  void registerWeightFunc(WeightFunc);
385  void registerErrorHandlerFunc(ErrorHandlerFunc);
386 
387 
388  // This is the library-side part of the CPodes constructor. This must
389  // be done prior to the client side construction.
390  void librarySideCPodesConstructor(ODEType, LinearMultistepMethod, NonlinearSystemIterationType);
391 
392  // Note that this routine MUST be called from client-side code so that
393  // it picks up exactly the static routines above which will agree with
394  // the client about the layout of the CPodesSystem virtual function table.
395  void clientSideCPodesConstructor() {
396  registerExplicitODEFunc(explicitODE_static);
397  registerImplicitODEFunc(implicitODE_static);
398  registerConstraintFunc(constraint_static);
399  registerProjectFunc(project_static);
400  registerQuadratureFunc(quadrature_static);
401  registerRootFunc(root_static);
402  registerWeightFunc(weight_static);
403  registerErrorHandlerFunc(errorHandler_static);
404  }
405 
406  // FOR INTERNAL USE ONLY
407 private:
408  class CPodesRep* rep;
409  friend class CPodesRep;
410 
411  const CPodesRep& getRep() const {assert(rep); return *rep;}
412  CPodesRep& updRep() {assert(rep); return *rep;}
413 
414  // Suppress copy constructor and default assigment operator.
415  CPodes(const CPodes&);
416  CPodes& operator=(const CPodes&);
417 };
418 
419 } // namespace SimTK
420 
421 #endif // SimTK_CPODES_H_
ODEType
Definition: SimTKcpodes.h:127
static int explicitODE_static(const CPodesSystem &sys, Real t, const Vector &y, Vector &fout)
Definition: SimTKcpodes.h:80
Definition: SimTKcpodes.h:176
that it contains events that have triggered but have not been processed Thus it is not a legitimate point along the systemТs trajectory and should never be returned to the caller For discussion we label this УimproperФ state EMBED Equation DSMT4 The event handler will correct the state creating a modified state with time still at thigh but not triggering any events We label the modified state EMBED Equation DSMT4 Once the systemТs handler the state EMBED Equation DSMT4 can be output as part of the trajectory and used as the initial condition for the next continuous interval Thus the time stepper generates the sequence of legitimate trajectory points just prior to event immediately followed by EMBED Equation DSMT4 which is the time at which the but after the event handler has modified the state to deal with those events Consecutive intervals I and I will consist of trajectory points EMBED Equation DSMT4 with round brackets indicating that EMBED Equation DSMT4 is not part of the trajectory Event handlers Event handlers are solvers which are able to take a state in an resolve the set up appropriate event triggers for the next and report back to the time stepper the degree to which the system continuity has been so that the integrator can be reinitialized appropriately An event handler can also indicate that the simulation should be in which case the time stepper will return the final state to its caller and disallow further time stepping Other event types Not all events have to be localized There are several special case clock user interrupt control is returned to the time stepper The time stepper can then declare that a scheduled event has call the systemТs event and reinitialize the integrator if continuity has been violated Time advanced events occur whenever the integrator has advanced time that at the end of every successful internal integration step These are generally restricted to discrete variable updates which do not affect the continuous such as min max values used only for reporting Normally the integrator does not return control at the end of a step
Definition: SimmathUserGuide.doc:227
static int weight_static(const CPodesSystem &sys, const Vector &y, Vector &weights)
Definition: SimTKcpodes.h:105
This abstract class defines the system to be integrated with SimTK CPodes.
Definition: SimTKcpodes.h:49
virtual int quadrature(Real t, const Vector &y, Vector &qout) const
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition: Assembler.h:37
Definition: SimTKcpodes.h:135
static int project_static(const CPodesSystem &sys, Real t, const Vector &ycur, Vector &corr, Real epsProj, Vector &err)
Definition: SimTKcpodes.h:92
Definition: SimTKcpodes.h:167
static int quadrature_static(const CPodesSystem &sys, Real t, const Vector &y, Vector &qout)
Definition: SimTKcpodes.h:97
LinearMultistepMethod
Definition: SimTKcpodes.h:133
ProjectionNorm
Definition: SimTKcpodes.h:152
static int constraint_static(const CPodesSystem &sys, Real t, const Vector &y, Vector &cout)
Definition: SimTKcpodes.h:88
virtual int implicitODE(Real t, const Vector &y, const Vector &yp, Vector &fout) const
Definition: SimTKcpodes.h:160
SimTK_Real Real
This is the default compiled-in floating point type for SimTK, either float or double.
Definition: SimTKcommon/include/SimTKcommon/internal/common.h:593
Definition: SimTKcpodes.h:141
╨╧ рб▒ с ■  ╖ ╣ ■    │ ┤ ╡ ╢                                                                                                                                                                                                                                                                                                                                                                                                                                     ье┴ А ° ┐ ч bjbjcTcT ┌┘ │ ├ ╗ t          ╖ Я ┴ K K K D      П П П А Л2 Ф П Z╞ j J a n u a r y
Definition: Simmatrix.doc:5
ConstraintLinearity
Definition: SimTKcpodes.h:158
Definition: SimTKcpodes.h:147
Includes internal headers providing declarations for the basic SimTK Core classes, including Simmatrix.
virtual int explicitODE(Real t, const Vector &y, Vector &fout) const
virtual int weight(const Vector &y, Vector &weights) const
virtual ~CPodesSystem()
Definition: SimTKcpodes.h:51
This is a straightforward translation of the Sundials CPODES C interface into C++.
Definition: SimTKcpodes.h:122
virtual int root(Real t, const Vector &y, const Vector &yp, Vector &gout) const
StepMode
Definition: SimTKcpodes.h:172
NonlinearSystemIterationType
Definition: SimTKcpodes.h:139
Definition: SimTKcpodes.h:148
This is the header file that every Simmath compilation unit should include first. ...
Definition: SimTKcpodes.h:129
Definition: SimTKcpodes.h:168
Definition: SimTKcpodes.h:174
virtual int project(Real t, const Vector &ycur, Vector &corr, Real epsProj, Vector &err) const
Definition: SimTKcpodes.h:154
ToleranceType
Definition: SimTKcpodes.h:145
Definition: SimTKcpodes.h:175
static int implicitODE_static(const CPodesSystem &sys, Real t, const Vector &y, const Vector &yp, Vector &fout)
Definition: SimTKcpodes.h:84
ProjectionFactorizationType
Definition: SimTKcpodes.h:164
#define SimTK_SIMMATH_EXPORT
Definition: SimTKmath/include/simmath/internal/common.h:64
static int root_static(const CPodesSystem &sys, Real t, const Vector &y, const Vector &yp, Vector &gout)
Definition: SimTKcpodes.h:101
virtual void errorHandler(int error_code, const char *module, const char *function, char *msg) const
static void errorHandler_static(const CPodesSystem &sys, int error_code, const char *module, const char *function, char *msg)
Definition: SimTKcpodes.h:109
virtual int constraint(Real t, const Vector &y, Vector &cout) const
Definition: SimTKcpodes.h:166