Simbody 3.7
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
38namespace SimTK {
39
50public:
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.
80static int explicitODE_static(const CPodesSystem& sys,
81 Real t, const Vector& y, Vector& fout)
82 { return sys.explicitODE(t,y,fout); }
83
84static 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
88static int constraint_static(const CPodesSystem& sys,
89 Real t, const Vector& y, Vector& cout)
90 { return sys.constraint(t,y,cout); }
91
92static 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
97static int quadrature_static(const CPodesSystem& sys,
98 Real t, const Vector& y, Vector& qout)
99 { return sys.quadrature(t,y,qout); }
100
101static 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
105static int weight_static(const CPodesSystem& sys,
106 const Vector& y, Vector& weights)
107 { return sys.weight(y,weights); }
108
109static 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
123public:
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
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.
217 Real t0, const Vector& y0, const Vector& yp0,
218 ToleranceType tt, Real reltol, void* abstol);
219
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.
227 const Vector& ctol);
228
229 // This tells CPodes to make use of the user's project()
230 // method from CPodesSystem.
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.
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) ;
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);
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);
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
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);
305 int getTolScaleFactor(Real* tolsfac);
306 int getErrWeights(Vector& eweight);
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);
344
345private:
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
407private:
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 assignment operator.
415 CPodes(const CPodes&);
416 CPodes& operator=(const CPodes&);
417};
418
419} // namespace SimTK
420
421#endif // SimTK_CPODES_H_
Includes internal headers providing declarations for the basic SimTK Core classes,...
This is the header file that every Simmath compilation unit should include first.
#define SimTK_SIMMATH_EXPORT
Definition: SimTKmath/include/simmath/internal/common.h:64
This abstract class defines the system to be integrated with SimTK CPodes.
Definition: SimTKcpodes.h:49
virtual int root(Real t, const Vector &y, const Vector &yp, Vector &gout) const
virtual int weight(const Vector &y, Vector &weights) const
virtual void errorHandler(int error_code, const char *module, const char *function, char *msg) const
virtual int explicitODE(Real t, const Vector &y, Vector &fout) const
virtual int constraint(Real t, const Vector &y, Vector &cout) const
virtual int project(Real t, const Vector &ycur, Vector &corr, Real epsProj, Vector &err) const
virtual ~CPodesSystem()
Definition: SimTKcpodes.h:51
virtual int implicitODE(Real t, const Vector &y, const Vector &yp, Vector &fout) const
virtual int quadrature(Real t, const Vector &y, Vector &qout) const
This is a straightforward translation of the Sundials CPODES C interface into C++.
Definition: SimTKcpodes.h:122
int getProjNumFailures(int *nprf)
int getNumFctEvals(int *nfevals)
int lapackBand(int N, int mupper, int mlower)
int getIntegratorStats(int *nsteps, int *nfevals, int *nlinsetups, int *netfails, int *qlast, int *qcur, Real *hinused, Real *hlast, Real *hcur, Real *tcur)
int lapackDense(int N)
int quadInit(const Vector &q0)
NonlinearSystemIterationType
Definition: SimTKcpodes.h:139
@ Newton
Definition: SimTKcpodes.h:141
int getRootInfo(int *rootsfound)
int reInit(CPodesSystem &sys, Real t0, const Vector &y0, const Vector &yp0, ToleranceType tt, Real reltol, void *abstol)
int setTolerances(int tol_type, Real reltol, void *abstol)
int getNumGEvals(int *ngevals)
int getNumNonlinSolvConvFails(int *nncfails)
ProjectionFactorizationType
Definition: SimTKcpodes.h:164
@ ProjectWithQR
Definition: SimTKcpodes.h:167
@ ProjectWithSchurComplement
Definition: SimTKcpodes.h:168
@ ProjectWithLU
Definition: SimTKcpodes.h:166
int getCurrentTime(Real *tcur)
LinearMultistepMethod
Definition: SimTKcpodes.h:133
@ BDF
Definition: SimTKcpodes.h:135
int getNonlinSolvStats(int *nniters, int *nncfails)
int setMinStep(Real hmin)
int dlsProjSetJacFn(void *jacP, void *jacP_data)
int dlsProjGetNumJacEvals(int *njPevals)
int setRootDirection(Array_< int > &rootdir)
ConstraintLinearity
Definition: SimTKcpodes.h:158
@ Linear
Definition: SimTKcpodes.h:160
int dlsSetJacFn(void *jac, void *jac_data)
int getNumLinSolvSetups(int *nlinsetups)
int dlsProjGetNumFctEvals(int *ncevalsLS)
int setStopTime(Real tstop)
int getCurrentOrder(int *qcur)
int getQuadNumFunEvals(int *nqevals)
int step(Real tout, Real *tret, Vector &y_inout, Vector &yp_inout, StepMode=Normal)
int getQuad(Real t, Vector &yQout)
int getWorkSpace(int *lenrw, int *leniw)
int setNonlinConvCoef(Real nlscoef)
int getNumNonlinSolvIters(int *nniters)
ToleranceType
Definition: SimTKcpodes.h:145
@ ScalarScalar
Definition: SimTKcpodes.h:147
@ ScalarVector
Definition: SimTKcpodes.h:148
int dlsGetWorkSpace(int *lenrwLS, int *leniwLS)
int setMaxStep(Real hmax)
int getProjNumCnstrEvals(int *nce)
ODEType
Definition: SimTKcpodes.h:127
@ ExplicitODE
Definition: SimTKcpodes.h:129
char * getReturnFlagName(int flag)
int setProjUpdateErrEst(bool proj_err)
int getDky(Real t, int k, Vector &dky)
int getNumSteps(int *nsteps)
int setProjFrequency(int proj_freq)
int getNumErrTestFails(int *netfails)
StepMode
Definition: SimTKcpodes.h:172
@ NormalTstop
Definition: SimTKcpodes.h:176
@ Normal
Definition: SimTKcpodes.h:174
@ OneStep
Definition: SimTKcpodes.h:175
int getNumStabLimOrderReds(int *nslred)
int setProjNonlinConvCoef(Real prjcoef)
int setProjLsetupFreq(int proj_lset_freq)
int getProjStats(int *nproj, int *nce, int *nsetupsP, int *nprf)
int projInit(ProjectionNorm, ConstraintLinearity, const Vector &ctol)
int getEstLocalErrors(Vector &ele)
int init(CPodesSystem &sys, Real t0, const Vector &y0, const Vector &yp0, ToleranceType tt, Real reltol, void *abstol)
int dlsGetLastFlag(int *flag)
ProjectionNorm
Definition: SimTKcpodes.h:152
@ L2Norm
Definition: SimTKcpodes.h:154
int dlsGetNumJacEvals(int *njevals)
int setMaxHnilWarns(int mxhnil)
int lapackDenseProj(int Nc, int Ny, ProjectionFactorizationType)
CPodes(ODEType ode=UnspecifiedODEType, LinearMultistepMethod lmm=UnspecifiedLinearMultistepMethod, NonlinearSystemIterationType nls=UnspecifiedNonlinearSystemIterationType)
Definition: SimTKcpodes.h:181
int setProjTestCnstr(bool test_cnstr)
int setMaxNumSteps(int mxsteps)
int setMaxNonlinIters(int maxcor)
int getLastOrder(int *qlast)
int setErrFile(FILE *errfp)
int setInitStep(Real hin)
int getProjNumProj(int *nproj)
int setQuadErrCon(bool errconQ, int tol_typeQ, Real reltolQ, void *abstolQ)
int getQuadErrWeights(Vector &eQweight)
int setMaxErrTestFails(int maxnef)
int getCurrentStep(Real *hcur)
int getRootWindow(Real *tLo, Real *tHi)
int setMaxConvFails(int maxncf)
int setMaxOrd(int maxord)
int setErrHandlerFn()
int getProjNumLinSolvSetups(int *nsetupsP)
int quadReInit(const Vector &q0)
int dlsGetNumFctEvals(int *nfevalsLS)
int setStabLimDet(bool stldet)
int getErrWeights(Vector &eweight)
int rootInit(int nrtfn)
char * dlsGetReturnFlagName(int flag)
int getQuadDky(Real t, int k, Vector &dky)
int getTolScaleFactor(Real *tolsfac)
int getLastStep(Real *hlast)
int getActualInitStep(Real *hinused)
This is the top-level SimTK namespace into which all SimTK names are placed to avoid collision with o...
Definition: Assembler.h:37
static int weight_static(const CPodesSystem &sys, const Vector &y, Vector &weights)
Definition: SimTKcpodes.h:105
static int explicitODE_static(const CPodesSystem &sys, Real t, const Vector &y, Vector &fout)
Definition: SimTKcpodes.h:80
static int quadrature_static(const CPodesSystem &sys, Real t, const Vector &y, Vector &qout)
Definition: SimTKcpodes.h:97
static int root_static(const CPodesSystem &sys, Real t, const Vector &y, const Vector &yp, Vector &gout)
Definition: SimTKcpodes.h:101
static int project_static(const CPodesSystem &sys, Real t, const Vector &ycur, Vector &corr, Real epsProj, Vector &err)
Definition: SimTKcpodes.h:92
static int constraint_static(const CPodesSystem &sys, Real t, const Vector &y, Vector &cout)
Definition: SimTKcpodes.h:88
static void errorHandler_static(const CPodesSystem &sys, int error_code, const char *module, const char *function, char *msg)
Definition: SimTKcpodes.h:109
static int implicitODE_static(const CPodesSystem &sys, Real t, const Vector &y, const Vector &yp, Vector &fout)
Definition: SimTKcpodes.h:84
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:606