Chaste  Release::3.4
CvodeAdaptor.cpp
1 /*
2 
3 Copyright (c) 2005-2016, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #ifdef CHASTE_CVODE
37 
38 #include "CvodeAdaptor.hpp"
39 
40 #include "Exception.hpp"
41 #include "TimeStepper.hpp"
43 #include "MathsCustomFunctions.hpp" // For tolerance check.
44 
45 #include <iostream>
46 #include <sstream>
47 
48 // CVODE headers
49 #include <sundials/sundials_nvector.h>
50 #include <cvode/cvode_dense.h>
51 
52 
69 int CvodeRhsAdaptor(realtype t, N_Vector y, N_Vector ydot, void* pData)
70 {
71  assert(pData != NULL);
72  CvodeData* p_data = (CvodeData*) pData;
73  // Get y, ydot into std::vector<>s
74  static std::vector<realtype> ydot_vec;
75  CopyToStdVector(y, *p_data->pY);
76  CopyToStdVector(ydot, ydot_vec);
77  // Call our function
78  try
79  {
80  p_data->pSystem->EvaluateYDerivatives(t, *(p_data->pY), ydot_vec);
81  }
82  catch (const Exception &e)
83  {
84  std::cerr << "CVODE RHS Exception: " << e.GetMessage() << std::endl << std::flush;
85  return -1;
86  }
87  // Copy derivative back
88  CopyFromStdVector(ydot_vec, ydot);
89  return 0;
90 }
91 
112 int CvodeRootAdaptor(realtype t, N_Vector y, realtype* pGOut, void* pData)
113 {
114  assert(pData != NULL);
115  CvodeData* p_data = (CvodeData*) pData;
116  // Get y into a std::vector
117  CopyToStdVector(y, *p_data->pY);
118  // Call our function
119  try
120  {
121  *pGOut = p_data->pSystem->CalculateRootFunction(t, *p_data->pY);
122  }
123  catch (const Exception &e)
124  {
125  std::cerr << "CVODE Root Exception: " << e.GetMessage() << std::endl << std::flush;
126  return -1;
127  }
128  return 0;
129 }
130 
131 // /**
132 // * Jacobian computation adaptor function.
133 // *
134 // * If solving an AbstractOdeSystemWithAnalyticJacobian, this function
135 // * can be used to allow CVODE to compute the Jacobian analytically.
136 // *
137 // * Note to self: can test using pSystem->GetUseAnalyticJacobian().
138 // */
139 // int CvodeDenseJacobianAdaptor(long int numberOfStateVariables, DenseMat J,
140 // realtype t, N_Vector y, N_Vector fy,
141 // void* pData,
142 // N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
143 // {
144 // AbstractOdeSystemWithAnalyticJacobian* pSystem
145 // = (AbstractOdeSystemWithAnalyticJacobian*) pData;
146 // // Get the current time step
147 // double dt;
148 // CVodeGetCurrentStep(CvodeMem, &dt);
149 // // Get std::vector<> for y and double** for J
150 // std::vector<realtype>& y_vec = *NV_DATA_STL(y);
151 // double** ppJ = J->data; // organised column-wise: J_{i,j} = ppJ[j][i]
152 // // Call our function
153 // try
154 // {
155 // pSystem->AnalyticJacobian(y_vec, ppJ, t, dt);
156 // }
157 // catch (const Exception &e)
158 // {
159 // std::cerr << "CVODE J Exception: " << e.GetMessage() << std::endl << std::flush;
160 // return -1;
161 // }
162 // // Update J (if needed)
163 // return 0;
164 // }
165 
166 
167 void CvodeErrorHandler(int errorCode, const char *module, const char *function,
168  char *message, void* pData)
169 {
170  std::stringstream err;
171  err << "CVODE Error " << errorCode << " in module " << module
172  << " function " << function << ": " << message;
173  std::cerr << "*" << err.str() << std::endl << std::flush;
174  // Throwing the exception here causes termination on Maverick (g++ 4.4)
175  //EXCEPTION(err.str());
176 }
177 
178 
180  std::vector<double>& rInitialY,
181  double startTime,
182  double maxStep)
183 {
184  assert(rInitialY.size() == pOdeSystem->GetNumberOfStateVariables());
185  assert(maxStep > 0.0);
186 
188  N_Vector initial_values = N_VMake_Serial(rInitialY.size(), &(rInitialY[0]));
189  assert(NV_DATA_S(initial_values) == &(rInitialY[0]));
190  assert(!NV_OWN_DATA_S(initial_values));
191 // std::cout << " Initial values: "; N_VPrint_Serial(initial_values);
192 // std::cout << " Rtol: " << mRelTol << ", Atol: " << mAbsTol << std::endl;
193 // std::cout << " Start: " << startTime << " max dt=" << maxStep << std::endl << std::flush;
194 
195  // Find out if we need to (re-)initialise
197  if (!reinit && !mForceMinimalReset)
198  {
199  const unsigned size = GetVectorSize(rInitialY);
200  for (unsigned i=0; i<size; i++)
201  {
203  {
204  reinit = true;
205  break;
206  }
207  }
208  }
209 
210  if (!mpCvodeMem) // First run of this solver, set up CVODE memory
211  {
212  // Set up CVODE's memory.
213  mpCvodeMem = CVodeCreate(CV_BDF, CV_NEWTON);
214  if (mpCvodeMem == NULL) EXCEPTION("Failed to SetupCvode CVODE"); // In one line to avoid coverage problem!
215  // Set error handler
216  CVodeSetErrHandlerFn(mpCvodeMem, CvodeErrorHandler, NULL);
217  // Set the user data
218  mData.pSystem = pOdeSystem;
219  mData.pY = &rInitialY;
220  #if CHASTE_SUNDIALS_VERSION >= 20400
221  CVodeSetUserData(mpCvodeMem, (void*)(&mData));
222  #else
223  CVodeSetFdata(mpCvodeMem, (void*)(&mData));
224  #endif
225 
226  // Setup CVODE
227  #if CHASTE_SUNDIALS_VERSION >= 20400
228  CVodeInit(mpCvodeMem, CvodeRhsAdaptor, startTime, initial_values);
229  CVodeSStolerances(mpCvodeMem, mRelTol, mAbsTol);
230  #else
231  CVodeMalloc(mpCvodeMem, CvodeRhsAdaptor, startTime, initial_values,
232  CV_SS, mRelTol, &mAbsTol);
233  #endif
234 
235  // Set the rootfinder function if wanted
236  if (mCheckForRoots)
237  {
238  #if CHASTE_SUNDIALS_VERSION >= 20400
239  CVodeRootInit(mpCvodeMem, 1, CvodeRootAdaptor);
240  #else
241  CVodeRootInit(mpCvodeMem, 1, CvodeRootAdaptor, (void*)(&mData));
242  #endif
243  }
244  // Attach a linear solver for Newton iteration
245  CVDense(mpCvodeMem, rInitialY.size());
246  }
247  else if (reinit) // Could be new ODE system, or new Y values
248  {
249  // Set the user data
250  mData.pSystem = pOdeSystem; // stays the same on a re-initialize
251  mData.pY = &rInitialY; // changes on a re-initialize
252  #if CHASTE_SUNDIALS_VERSION >= 20400
253  CVodeSetUserData(mpCvodeMem, (void*)(&mData));
254  #else
255  CVodeSetFdata(mpCvodeMem, (void*)(&mData));
256  #endif
257 
258  #if CHASTE_SUNDIALS_VERSION >= 20400
259  CVodeReInit(mpCvodeMem, startTime, initial_values);
260  CVodeSStolerances(mpCvodeMem, mRelTol, mAbsTol);
261  #else
262  CVodeReInit(mpCvodeMem, CvodeRhsAdaptor, startTime, initial_values,
263  CV_SS, mRelTol, &mAbsTol);
264  #endif
265 
266  // Attach a linear solver for Newton iteration
267  CVDense(mpCvodeMem, rInitialY.size());
268  }
269 
270  CVodeSetMaxStep(mpCvodeMem, maxStep);
271  // Change max steps if wanted
272  if (mMaxSteps > 0)
273  {
274  CVodeSetMaxNumSteps(mpCvodeMem, mMaxSteps);
275  CVodeSetMaxErrTestFails(mpCvodeMem, 15);
276  }
277  DeleteVector(initial_values);
278 }
279 
281 {
282  if (mpCvodeMem)
283  {
284  CVodeFree(&mpCvodeMem);
285  }
286  mpCvodeMem = NULL;
287 }
288 
289 
290 void CvodeAdaptor::SetForceReset(bool autoReset)
291 {
292  mForceReset = autoReset;
293  if (mForceReset)
294  {
295  SetMinimalReset(false);
296  ResetSolver();
297  }
298 }
299 
300 void CvodeAdaptor::SetMinimalReset(bool minimalReset)
301 {
302  mForceMinimalReset = minimalReset;
303  if (mForceMinimalReset)
304  {
305  SetForceReset(false);
306  }
307 }
308 
310 {
311  // Doing this makes the SetupCvode think it needs to reset things.
313 }
314 
315 void CvodeAdaptor::CvodeError(int flag, const char * msg)
316 {
317  std::stringstream err;
318  char* p_flag_name = CVodeGetReturnFlagName(flag);
319  err << msg << ": " << p_flag_name;
320  free(p_flag_name);
321  std::cerr << err.str() << std::endl << std::flush;
322  EXCEPTION(err.str());
323 }
324 
325 
327  std::vector<double>& rYValues,
328  double startTime,
329  double endTime,
330  double maxStep,
331  double timeSampling)
332 {
333  assert(endTime > startTime);
334  assert(timeSampling > 0.0);
335 
336  mStoppingEventOccurred = false;
337  if (mCheckForRoots && pOdeSystem->CalculateStoppingEvent(startTime, rYValues) == true)
338  {
339  EXCEPTION("(Solve with sampling) Stopping event is true for initial condition");
340  }
341 
342  SetupCvode(pOdeSystem, rYValues, startTime, maxStep);
343 
344  TimeStepper stepper(startTime, endTime, timeSampling);
345  N_Vector yout = N_VMake_Serial(rYValues.size(), &(rYValues[0]));
346 
347  // Set up ODE solution
348  OdeSolution solutions;
349  solutions.SetNumberOfTimeSteps(stepper.EstimateTimeSteps());
350  solutions.rGetSolutions().push_back(rYValues);
351  solutions.rGetTimes().push_back(startTime);
352  solutions.SetOdeSystemInformation(pOdeSystem->GetSystemInformation());
353 
354  // Main time sampling loop
355  while (!stepper.IsTimeAtEnd() && !mStoppingEventOccurred)
356  {
357  // This should stop CVODE going past the end of where we wanted and interpolating back.
358  int ierr = CVodeSetStopTime(mpCvodeMem, stepper.GetNextTime());
359  assert(ierr == CV_SUCCESS); UNUSED_OPT(ierr); // avoid unused var warning
360 
361  double tend;
362  ierr = CVode(mpCvodeMem, stepper.GetNextTime(), yout, &tend, CV_NORMAL);
363  if (ierr<0)
364  {
365  FreeCvodeMemory();
366  DeleteVector(yout);
367  CvodeError(ierr, "CVODE failed to solve system");
368  }
369  // Store solution
370  solutions.rGetSolutions().push_back(rYValues);
371  solutions.rGetTimes().push_back(tend);
372  if (ierr == CV_ROOT_RETURN)
373  {
374  // Stopping event occurred
375  mStoppingEventOccurred = true;
376  mStoppingTime = tend;
377  }
378  mLastSolutionTime = tend;
379  stepper.AdvanceOneTimeStep();
380  }
381 
382  // stepper.EstimateTimeSteps may have been an overestimate...
383  solutions.SetNumberOfTimeSteps(stepper.GetTotalTimeStepsTaken());
384 
385  int ierr = CVodeGetLastStep(mpCvodeMem, &mLastInternalStepSize);
386  assert(ierr == CV_SUCCESS); UNUSED_OPT(ierr); // avoid unused var warning
388  DeleteVector(yout);
389 
390  return solutions;
391 }
392 
393 
395  std::vector<double>& rYValues,
396  double startTime,
397  double endTime,
398  double maxStep)
399 {
400  assert(endTime > startTime);
401 
402  mStoppingEventOccurred = false;
403  if (mCheckForRoots && pOdeSystem->CalculateStoppingEvent(startTime, rYValues) == true)
404  {
405  EXCEPTION("(Solve) Stopping event is true for initial condition");
406  }
407 
408  SetupCvode(pOdeSystem, rYValues, startTime, maxStep);
409 
410  N_Vector yout = N_VMake_Serial(rYValues.size(), &(rYValues[0]));
411 
412  // This should stop CVODE going past the end of where we wanted and interpolating back.
413  int ierr = CVodeSetStopTime(mpCvodeMem, endTime);
414  assert(ierr == CV_SUCCESS); UNUSED_OPT(ierr); // avoid unused var warning
415 
416  double tend;
417  ierr = CVode(mpCvodeMem, endTime, yout, &tend, CV_NORMAL);
418  if (ierr<0)
419  {
420  FreeCvodeMemory();
421  DeleteVector(yout);
422  CvodeError(ierr, "CVODE failed to solve system");
423  }
424  if (ierr == CV_ROOT_RETURN)
425  {
426  // Stopping event occurred
427  mStoppingEventOccurred = true;
428  mStoppingTime = tend;
429  }
430  assert(NV_DATA_S(yout) == &(rYValues[0]));
431  assert(!NV_OWN_DATA_S(yout));
432 
433 // long int steps;
434 // CVodeGetNumSteps(mpCvodeMem, &steps);
435 // std::cout << " Solved to " << endTime << " in " << steps << " steps.\n";
436 
437  ierr = CVodeGetLastStep(mpCvodeMem, &mLastInternalStepSize);
438  assert(ierr == CV_SUCCESS); UNUSED_OPT(ierr); // avoid unused var warning
439  RecordStoppingPoint(tend, yout);
440  DeleteVector(yout);
441 }
442 
443 CvodeAdaptor::CvodeAdaptor(double relTol, double absTol)
445  mpCvodeMem(NULL),
446  mRelTol(relTol),
447  mAbsTol(absTol),
448  mLastInternalStepSize(-0.0),
449  mMaxSteps(0),
450  mCheckForRoots(false),
451  mLastSolutionState(NULL),
452  mLastSolutionTime(0.0),
453 #if CHASTE_SUNDIALS_VERSION >= 20400
454  mForceReset(false),
455 #else
456  mForceReset(true),
457 #endif
458  mForceMinimalReset(false)
459 {
460 }
461 
463 {
464  FreeCvodeMemory();
466 }
467 
468 void CvodeAdaptor::RecordStoppingPoint(double stopTime, N_Vector yEnd)
469 {
470  if (!mForceReset)
471  {
472  const unsigned size = GetVectorSize(yEnd);
474  for (unsigned i=0; i<size; i++)
475  {
477  }
478  mLastSolutionTime = stopTime;
479  }
480 }
481 
482 void CvodeAdaptor::SetTolerances(double relTol, double absTol)
483 {
484  mRelTol = relTol;
485  mAbsTol = absTol;
486 }
487 
489 {
490  return mRelTol;
491 }
492 
494 {
495  return mAbsTol;
496 }
497 
499 {
500  return mLastInternalStepSize;
501 }
502 
504 {
505  mCheckForRoots = true;
506 }
507 
508 void CvodeAdaptor::SetMaxSteps(long int numSteps)
509 {
510  mMaxSteps = numSteps;
511 }
512 
514 {
515  return mMaxSteps;
516 }
517 
518 // Serialization for Boost >= 1.36
521 
522 #endif // CHASTE_CVODE
AbstractOdeSystem * pSystem
void SetupCvode(AbstractOdeSystem *pOdeSystem, std::vector< double > &rInitialY, double startTime, double maxStep)
static bool WithinAnyTolerance(double number1, double number2, double relTol=DBL_EPSILON, double absTol=DBL_EPSILON, bool printError=false)
CvodeAdaptor(double relTol=1e-4, double absTol=1e-6)
std::vector< std::vector< double > > & rGetSolutions()
void ResetSolver()
void CheckForStoppingEvents()
#define EXCEPTION(message)
Definition: Exception.hpp:143
N_Vector mLastSolutionState
virtual void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)=0
void FreeCvodeMemory()
long int GetMaxSteps()
double mLastInternalStepSize
void AdvanceOneTimeStep()
bool IsTimeAtEnd() const
double GetAbsoluteTolerance()
void CreateVectorIfEmpty(VECTOR &rVec, unsigned size)
void SetOdeSystemInformation(boost::shared_ptr< const AbstractOdeSystemInformation > pOdeSystemInfo)
Definition: OdeSolution.cpp:66
void CvodeError(int flag, const char *msg)
void SetMinimalReset(bool minimalReset)
std::string GetMessage() const
Definition: Exception.cpp:87
bool mForceMinimalReset
virtual double CalculateRootFunction(double time, const std::vector< double > &rY)
void SetVectorComponent(VECTOR &rVec, unsigned index, double value)
std::vector< realtype > * pY
unsigned GetTotalTimeStepsTaken() const
boost::shared_ptr< const AbstractOdeSystemInformation > GetSystemInformation() const
std::vector< double > & rGetTimes()
unsigned EstimateTimeSteps() const
void SetTolerances(double relTol=1e-4, double absTol=1e-6)
#define UNUSED_OPT(var)
Definition: Exception.hpp:216
#define CHASTE_CLASS_EXPORT(T)
void RecordStoppingPoint(double stopTime, N_Vector yEnd)
void SetForceReset(bool autoReset)
void CopyFromStdVector(const std::vector< double > &rSrc, VECTOR &rDest)
double mLastSolutionTime
CvodeData mData
void CopyToStdVector(const VECTOR &rSrc, std::vector< double > &rDest)
long int mMaxSteps
double GetNextTime() const
double GetVectorComponent(const VECTOR &rVec, unsigned index)
void SetMaxSteps(long int numSteps)
double GetLastStepSize()
double GetRelativeTolerance()
OdeSolution Solve(AbstractOdeSystem *pOdeSystem, std::vector< double > &rYValues, double startTime, double endTime, double maxStep, double timeSampling)
void SetNumberOfTimeSteps(unsigned numTimeSteps)
Definition: OdeSolution.cpp:58
void DeleteVector(VECTOR &rVec)
unsigned GetVectorSize(const VECTOR &rVec)
virtual bool CalculateStoppingEvent(double time, const std::vector< double > &rY)