Chaste  Release::3.4
TimeStepper.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 #include <cassert>
37 #include <cmath>
38 
39 #include "TimeStepper.hpp"
40 #include "Exception.hpp"
41 #include "MathsCustomFunctions.hpp"
42 
43 TimeStepper::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  }
134  mTime = mNextTime;
135 
137 }
138 
139 double TimeStepper::GetTime() const
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 {
177  return mTotalTimeStepsTaken;
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 }
double GetNextTimeStep()
double GetTime() const
#define EXCEPTION(message)
Definition: Exception.hpp:143
bool Divides(double smallerNumber, double largerNumber)
void ResetTimeStep(double dt)
void AdvanceOneTimeStep()
bool IsTimeAtEnd() const
unsigned GetTotalTimeStepsTaken() const
double GetIdealTimeStep()
double CalculateNextTime()
unsigned EstimateTimeSteps() const
double mNextTime
std::vector< double > mAdditionalTimesDeprecated
double GetNextTime() const
unsigned mTotalTimeStepsTaken
double mEpsilon