Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
TimeStepper.cpp
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#include <cassert>
37#include <cmath>
38
39#include "TimeStepper.hpp"
40#include "Exception.hpp"
42
43TimeStepper::TimeStepper(double startTime, double endTime, double dt, bool enforceConstantTimeStep, std::vector<double> additionalTimes)
44 : mStart(startTime),
45 mEnd(endTime),
46 mDt(dt),
47 mTotalTimeStepsTaken(0),
48 mAdditionalTimesReachedDeprecated(0),
49 mTime(startTime),
50 mEpsilon(DBL_EPSILON)
51{
52 if (startTime > endTime)
53 {
54 EXCEPTION("The simulation duration must be positive, not " << endTime-startTime);
55 }
56
57 // Remove any additionalTimes entries which fall too close to a time when the stepper would stop anyway
58 for (unsigned i=0; i<additionalTimes.size(); i++)
59 {
60 if (i > 0)
61 {
62 if (additionalTimes[i-1] >= additionalTimes[i])
63 {
64 EXCEPTION("The additional times vector should be in ascending numerical order; "
65 "entry " << i << " is less than or equal to entry " << i-1 << ".");
66 }
67 }
68
69 double time_interval = additionalTimes[i] - startTime;
70
71 // When mDt divides this interval (and the interval is positive) then we are going there anyway
72 if (!Divides(mDt, time_interval) && (time_interval > DBL_EPSILON))
73 {
74 //mAdditionalTimes.push_back(additionalTimes[i]);
75 EXCEPTION("Additional times are now deprecated. Use only to check whether the given times are met: e.g. Electrode events should only happen on printing steps.");
76 }
77 }
78
79 /*
80 * Note that when mEnd is large then the error of subtracting two numbers of
81 * that magnitude is about DBL_EPSILON*mEnd (1e-16*mEnd). When mEnd is small
82 * then the error should be around DBL_EPSILON.
83 */
84 if (mEnd > 1.0)
85 {
86 mEpsilon = DBL_EPSILON*mEnd;
87 }
88
89 // If enforceConstantTimeStep check whether the times are such that we won't have a variable dt
90 if (enforceConstantTimeStep)
91 {
92 double expected_end_time = mStart + mDt*EstimateTimeSteps();
93
94 if (fabs( mEnd - expected_end_time ) > mEpsilon)
95 {
96 EXCEPTION("TimeStepper estimates non-constant timesteps will need to be used: check timestep "
97 "divides (end_time-start_time) (or divides printing timestep). "
98 "[End time=" << mEnd << "; start=" << mStart << "; dt=" << mDt << "; error="
99 << fabs(mEnd-expected_end_time) << "]");
100 }
101 }
102
104}
105
107{
108 double next_time = mStart + (mTotalTimeStepsTaken + 1)*mDt;
109
110 // Does the next time bring us very close to the end time?
111 // Note that the inequality in this guard matches the inversion of the guard in the enforceConstantTimeStep
112 // calculation of the constructor
113 if (mEnd - next_time <= mEpsilon)
114 {
115 next_time = mEnd;
116 }
117
118 assert(mAdditionalTimesDeprecated.empty());
119
120 return next_time;
121}
122
124{
126 if (mTotalTimeStepsTaken == 0)
127 {
128 EXCEPTION("Time step counter has overflowed.");
129 }
130 if (mTime >= mNextTime)
131 {
132 EXCEPTION("TimeStepper incremented beyond end time.");
133 }
135
137}
138
140{
141 return mTime;
142}
143
145{
146 return mNextTime;
147}
148
150{
151 double dt = mDt;
152
153 if (mNextTime >= mEnd)
154 {
155 dt = mEnd - mTime;
156 }
157 assert(mAdditionalTimesDeprecated.empty());
158 return dt;
159}
161{
162 return mDt;
163}
164
166{
167 return (mTime >= mEnd);
168}
169
171{
172 return (unsigned) floor((mEnd - mStart)/mDt + 0.5);
173}
174
176{
178}
179
181{
182 assert(dt > 0);
183 /*
184 * The error in subtracting two numbers of the same magnitude is about
185 * DBL_EPSILON times that magnitude (we use the sum of the two numbers
186 * here as a conservative estimate of their maximum). When both mDt and
187 * dt are small then the error should be around DBL_EPSILON.
188 */
189 double scale = DBL_EPSILON*(mDt + dt);
190 if (mDt + dt < 1.0)
191 {
192 scale = DBL_EPSILON;
193 }
194 if (fabs(mDt-dt) > scale)
195 {
196 mDt = dt;
197 mStart = mTime;
199
201 }
202}
#define EXCEPTION(message)
unsigned mTotalTimeStepsTaken
void ResetTimeStep(double dt)
double GetNextTimeStep()
bool IsTimeAtEnd() const
std::vector< double > mAdditionalTimesDeprecated
unsigned GetTotalTimeStepsTaken() const
double GetTime() const
double CalculateNextTime()
double GetIdealTimeStep()
void AdvanceOneTimeStep()
double GetNextTime() const
unsigned EstimateTimeSteps() const