Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
AbstractCvodeSystem.hpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, 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
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#ifdef CHASTE_CVODE
37#ifndef _ABSTRACTCVODESYSTEM_HPP_
38#define _ABSTRACTCVODESYSTEM_HPP_
39
40#include <algorithm>
41#include <string>
42#include <vector>
43
44// This is only needed to prevent compilation errors on PETSc 2.2/Boost 1.33.1 combo
46
47// Chaste includes
48#include "AbstractParameterisedSystem.hpp"
49#include "Exception.hpp"
50#include "OdeSolution.hpp"
52
53// Serialiazation
54#include <boost/serialization/split_member.hpp>
55#include <boost/serialization/vector.hpp>
57#include "ClassIsAbstract.hpp"
58
59// CVODE headers
60#include <nvector/nvector_serial.h>
61
62#if CHASTE_SUNDIALS_VERSION >= 30000
63#include <cvode/cvode_direct.h> /* access to CVDls interface */
64#include <sundials/sundials_types.h> /* defs. of realtype, sunindextype */
65#include <sunlinsol/sunlinsol_dense.h> /* access to dense SUNLinearSolver */
66#include <sunmatrix/sunmatrix_dense.h> /* access to dense SUNMatrix */
67#else
68#include <sundials/sundials_dense.h> /* definitions DlsMat DENSE_ELEM */
69#endif
70
71// CVODE changed their dense matrix type...
72#if CHASTE_SUNDIALS_VERSION >= 30000
73#define CHASTE_CVODE_DENSE_MATRIX SUNMatrix
74#elif CHASTE_SUNDIALS_VERSION >= 20400
75#define CHASTE_CVODE_DENSE_MATRIX DlsMat
76#else
77#define CHASTE_CVODE_DENSE_MATRIX DenseMat
78#endif
79
80// CVODE changed their way of referencing elements of a matrix. So we will copy their notation in examples.
81#if CHASTE_SUNDIALS_VERSION >= 30000
82#define IJth(A, i, j) SM_ELEMENT_D(A, i, j)
83#else
84#define IJth(A, i, j) DENSE_ELEM(A, i, j)
85#endif
86
142{
143private:
144 friend class TestAbstractCvodeSystem;
145
146 friend class boost::serialization::access;
153 template <class Archive>
154 void save(Archive& archive, const unsigned int version) const
155 {
156 // Despite the fact that 3 of these variables actually live in our base class,
157 // we still archive them here to maintain backwards compatibility,
158 // this doesn't hurt
160 archive& mUseAnalyticJacobian;
161
162 if (version >= 1u)
163 {
164 archive& mHasAnalyticJacobian;
165 }
166
167 // Convert from N_Vector to std::vector for serialization
168 const std::vector<double> state_vars = MakeStdVec(mStateVariables);
169 archive& state_vars;
170 const std::vector<double> params = MakeStdVec(mParameters);
171 archive& params;
172 archive& rGetParameterNames();
173
174 archive& mLastSolutionTime;
175 archive& mForceReset;
176 archive& mForceMinimalReset;
177 archive& mRelTol;
178 archive& mAbsTol;
179 archive& mMaxSteps;
180 archive& mLastInternalStepSize;
181
182 // We don't bother archiving CVODE's internal data, because it is missing then we'll just
183 // get a new solver being initialised after a save/load.
184
185 // This is always set up by subclass constructors, and is essentially
186 // 'static' data, so shouldn't go in the archive.
187 //archive &mpSystemInfo;
188 }
195 template <class Archive>
196 void load(Archive& archive, const unsigned int version)
197 {
199 archive& mUseAnalyticJacobian;
200
201 // This is pretty much what the code was saying before.
203 if (version >= 1u)
204 { // Overwrite if it has been archived though
205 archive& mHasAnalyticJacobian;
206 }
207
208 std::vector<double> state_vars;
209 archive& state_vars;
211
212 std::vector<double> parameters;
213 archive& parameters;
214
215 std::vector<std::string> param_names;
216 archive& param_names;
217 archive& mLastSolutionTime;
218 archive& mForceReset;
219 archive& mForceMinimalReset;
220 archive& mRelTol;
221 archive& mAbsTol;
222 archive& mMaxSteps;
223 archive& mLastInternalStepSize;
224
225 // We don't bother archiving CVODE's internal data, because it is missing then we'll just
226 // get a new solver being initialised after a save/load.
227
228 // Do some checking on the parameters
229 CheckParametersOnLoad(parameters, param_names);
230 }
231 BOOST_SERIALIZATION_SPLIT_MEMBER()
232
233
240 void SetupCvode(N_Vector initialConditions,
241 realtype tStart,
242 realtype maxDt);
243
248 void RecordStoppingPoint(double stopTime);
249
251 void FreeCvodeMemory();
252
263 void CvodeError(int flag, const char* msg, const double& rTime,
264 const double& rStartTime, const double& rEndTime);
265
268
271
274
277
278#if CHASTE_SUNDIALS_VERSION >= 30000
280 SUNMatrix mpSundialsDenseMatrix;
282 SUNLinearSolver mpSundialsLinearSolver;
283#endif
284
285protected:
288
291
293 double mRelTol;
294
296 double mAbsTol;
297
300
305 long int mMaxSteps;
306
309
314 void Init();
315
316public:
322 AbstractCvodeSystem(unsigned numberOfStateVariables);
323
327 virtual ~AbstractCvodeSystem();
328
336 virtual void EvaluateYDerivatives(realtype time,
337 const N_Vector y,
338 N_Vector ydot)
339 = 0;
340
356 virtual void EvaluateAnalyticJacobian(realtype time, N_Vector y, N_Vector ydot,
357 CHASTE_CVODE_DENSE_MATRIX jacobian,
358 N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
359 {
360 EXCEPTION("No analytic Jacobian has been defined for this system.");
361 }
362
373 void SetForceReset(bool autoReset);
374
383 void SetMinimalReset(bool minimalReset);
384
391 bool GetMinimalReset();
392
399 bool GetForceReset();
400
411 void ResetSolver();
412
432 OdeSolution Solve(realtype tStart,
433 realtype tEnd,
434 realtype maxDt,
435 realtype tSamp);
436
452 void Solve(realtype tStart,
453 realtype tEnd,
454 realtype maxDt);
455
462 void SetMaxSteps(long int numSteps);
463
468 long int GetMaxSteps();
469
477 void SetTolerances(double relTol = 1e-5, double absTol = 1e-7);
478
482 double GetRelativeTolerance();
483
487 double GetAbsoluteTolerance();
488
492 double GetLastStepSize();
493
494 /* NB This needs making into a doxygen comment if you bring the method back in.
495 *
496 * An alternative approach to stopping events; currently only useful with CVODE.
497 * CVODE can search for roots (zeros) of this function while solving the ODE system,
498 * and home in on them to find sign transitions to high precision.
499 *
500 * The default implementation here fakes a root function using CalculateStoppingEvent.
501 *
502 * @param time the current time
503 * @param rY the current values of the state variables
504 */
505 // virtual double CalculateRootFunction(double time, const std::vector<double>& rY);
506
510 bool GetUseAnalyticJacobian() const;
511
515 bool HasAnalyticJacobian() const;
516
523 void ForceUseOfNumericalJacobian(bool useNumericalJacobian = true);
524
525 // The following method may be useful to identify problems with the Analytic Jacobians, if anything goes wrong,
526 // but #1795 seems to have got these working OK, so commented out for now.
527
528 // /*
529 // * Compare the calculated analytic jacobian to a numerical approximation, and throw if it looks silly.
530 // *
531 // * @param time the current time
532 // * @param y the current state variables
533 // * @param jacobian the analytic jacobian matrix
534 // * @param tmp1 working memory of the correct size provided by CVODE for temporary calculations
535 // * @param tmp2 working memory of the correct size provided by CVODE for temporary calculations
536 // * @param tmp3 working memory of the correct size provided by CVODE for temporary calculations
537 // */
538 // void CheckAnalyticJacobian(realtype time, N_Vector y, N_Vector ydot,
539 // CHASTE_CVODE_DENSE_MATRIX jacobian,
540 // N_Vector tmp1, N_Vector tmp2, N_Vector tmp3);
541};
542
544BOOST_CLASS_VERSION(AbstractCvodeSystem, 1u)
545
546#endif //_ABSTRACTCVODESYSTEM_HPP_
547#endif // CHASTE_CVODE
#define CLASS_IS_ABSTRACT(T)
#define EXCEPTION(message)
std::vector< double > MakeStdVec(N_Vector v)
void CopyFromStdVector(const std::vector< double > &rSrc, VECTOR &rDest)
void SetTolerances(double relTol=1e-5, double absTol=1e-7)
bool GetMinimalReset()
Get whether we want to run with minimal reset or not (no reinitialisation of the solver if variables ...
void SetForceReset(bool autoReset)
void load(Archive &archive, const unsigned int version)
void SetupCvode(N_Vector initialConditions, realtype tStart, realtype maxDt)
void CvodeError(int flag, const char *msg, const double &rTime, const double &rStartTime, const double &rEndTime)
void RecordStoppingPoint(double stopTime)
void SetMaxSteps(long int numSteps)
OdeSolution Solve(realtype tStart, realtype tEnd, realtype maxDt, realtype tSamp)
void ForceUseOfNumericalJacobian(bool useNumericalJacobian=true)
virtual void EvaluateYDerivatives(realtype time, const N_Vector y, N_Vector ydot)=0
void save(Archive &archive, const unsigned int version) const
virtual void EvaluateAnalyticJacobian(realtype time, N_Vector y, N_Vector ydot, CHASTE_CVODE_DENSE_MATRIX jacobian, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
void SetMinimalReset(bool minimalReset)
bool GetForceReset()
Get whether we will force a solver reset on every call to Solve()
void CheckParametersOnLoad(const std::vector< double > &rParameters, const std::vector< std::string > &rParameterNames)
const std::vector< std::string > & rGetParameterNames() const