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
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
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