Chaste  Release::2018.1
AbstractBackwardEulerCardiacCell.hpp
1 /*
2 
3 Copyright (c) 2005-2018, 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 
37 #ifndef ABSTRACTBACKWARDEULERCARDIACCELL_HPP_
38 #define ABSTRACTBACKWARDEULERCARDIACCELL_HPP_
39 
40 #include <cassert>
41 #include <cmath>
42 
43 #include "ChasteSerialization.hpp"
44 #include <boost/serialization/base_object.hpp>
45 #include "ClassIsAbstract.hpp"
46 
47 #include "AbstractCardiacCell.hpp"
48 #include "Exception.hpp"
49 #include "PetscTools.hpp"
50 #include "TimeStepper.hpp"
51 
66 template<unsigned SIZE>
68 {
69  private:
78  template<class Archive>
79  void serialize(Archive & archive, const unsigned int version)
80  {
81  // This calls serialize on the base class.
82  archive & boost::serialization::base_object<AbstractCardiacCell>(*this);
83  }
84 public:
85 
102  unsigned numberOfStateVariables,
103  unsigned voltageIndex,
104  boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus);
105 
108 
116  virtual void ComputeResidual(double time, const double rCurrentGuess[SIZE], double rResidual[SIZE])=0;
117 
125  virtual void ComputeJacobian(double time, const double rCurrentGuess[SIZE], double rJacobian[SIZE][SIZE])=0;
126 
139  OdeSolution Compute(double tStart, double tEnd, double tSamp=0.0);
140 
150  void ComputeExceptVoltage(double tStart, double tEnd);
151 
159  void SolveAndUpdateState(double tStart, double tEnd);
160 
161 private:
162 // LCOV_EXCL_START
170  void EvaluateYDerivatives(double time, const std::vector<double> &rY, std::vector<double> &rDY)
171  {
173  }
174 // LCOV_EXCL_STOP
175 
176 protected:
185  virtual void ComputeOneStepExceptVoltage(double tStart)=0;
186 
194  virtual void UpdateTransmembranePotential(double time)=0;
195 };
196 
197 
201 /*
202  * NOTE: Explicit instantiation is not used for this class, because the SIZE
203  * template parameter could take arbitrary values.
204  */
205 
206 
207 template <unsigned SIZE>
209  unsigned numberOfStateVariables,
210  unsigned voltageIndex,
211  boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
212  : AbstractCardiacCell(boost::shared_ptr<AbstractIvpOdeSolver>(),
213  numberOfStateVariables,
214  voltageIndex,
215  pIntracellularStimulus)
216 {}
217 
218 template <unsigned SIZE>
220 {}
221 
222 template <unsigned SIZE>
223 OdeSolution AbstractBackwardEulerCardiacCell<SIZE>::Compute(double tStart, double tEnd, double tSamp)
224 {
225  // In this method, we iterate over timesteps, doing the following for each:
226  // - update V using a forward Euler step
227  // - call ComputeExceptVoltage(t) to update the remaining state variables
228  // using backward Euler
229 
230  // Check length of time interval
231  if (tSamp < mDt)
232  {
233  tSamp = mDt;
234  }
235  double _n_steps = (tEnd - tStart) / tSamp;
236  const unsigned n_steps = (unsigned) floor(_n_steps+0.5);
237  assert(fabs(tStart+n_steps*tSamp - tEnd) < 1e-12);
238  const unsigned n_small_steps = (unsigned) floor(tSamp/mDt+0.5);
239  assert(fabs(mDt*n_small_steps - tSamp) < 1e-12);
240 
241  // Initialise solution store
242  OdeSolution solutions;
243  solutions.SetNumberOfTimeSteps(n_steps);
244  solutions.rGetSolutions().push_back(rGetStateVariables());
245  solutions.rGetTimes().push_back(tStart);
246  solutions.SetOdeSystemInformation(this->mpSystemInfo);
247 
248  // Loop over time
249  double curr_time = tStart;
250  for (unsigned i=0; i<n_steps; i++)
251  {
252  for (unsigned j=0; j<n_small_steps; j++)
253  {
254  curr_time = tStart + i*tSamp + j*mDt;
255 
256  // Compute next value of V
257  UpdateTransmembranePotential(curr_time);
258 
259  // Compute other state variables
260  ComputeOneStepExceptVoltage(curr_time);
261 
262  // check gating variables are still in range
263  VerifyStateVariables();
264  }
265 
266  // Update solutions
267  solutions.rGetSolutions().push_back(rGetStateVariables());
268  solutions.rGetTimes().push_back(curr_time+mDt);
269  }
270 
271  return solutions;
272 }
273 
274 template <unsigned SIZE>
276 {
277  // This method iterates over timesteps, calling ComputeExceptVoltage(t) at
278  // each one, to update all state variables except for V, using backward Euler.
279 
280  // Check length of time interval
281  unsigned n_steps = (unsigned)((tEnd - tStart) / mDt + 0.5);
282  assert(fabs(tStart + n_steps*mDt - tEnd) < 1e-12);
283 
284  // Loop over time
285  double curr_time;
286  for (unsigned i=0; i<n_steps; i++)
287  {
288  curr_time = tStart + i*mDt;
289 
290  // Compute other state variables
291  ComputeOneStepExceptVoltage(curr_time);
292 
293 #ifndef NDEBUG
294  // Check gating variables are still in range
295  VerifyStateVariables();
296 #endif // NDEBUG
297  }
298 }
299 
300 template<unsigned SIZE>
302 {
303  TimeStepper stepper(tStart, tEnd, mDt);
304 
305  while (!stepper.IsTimeAtEnd())
306  {
307  double time = stepper.GetTime();
308 
309  // Compute next value of V
310  UpdateTransmembranePotential(time);
311 
312  // Compute other state variables
313  ComputeOneStepExceptVoltage(time);
314 
315  // check gating variables are still in range
316  VerifyStateVariables();
317 
318  stepper.AdvanceOneTimeStep();
319  }
320 }
321 
325 
329 template<>
331 {
332 private:
334  friend class boost::serialization::access;
341  template<class Archive>
342  void serialize(Archive & archive, const unsigned int version)
343  {
344  // This calls serialize on the base class.
345  archive & boost::serialization::base_object<AbstractCardiacCell>(*this);
346  }
347 
348 public:
357  AbstractBackwardEulerCardiacCell(unsigned numberOfStateVariables,
358  unsigned voltageIndex,
359  boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
360  : AbstractCardiacCell(boost::shared_ptr<AbstractIvpOdeSolver>(),
361  numberOfStateVariables,
362  voltageIndex,
363  pIntracellularStimulus)
364  {}
365 
368  {}
369 
382  OdeSolution Compute(double tStart, double tEnd, double tSamp=0.0)
383  {
384  // In this method, we iterate over timesteps, doing the following for each:
385  // - update V using a forward Euler step
386  // - call ComputeExceptVoltage(t) to update the remaining state variables
387  // using backward Euler
388 
389  // Check length of time interval
390  if (tSamp < mDt)
391  {
392  tSamp = mDt;
393  }
394  double _n_steps = (tEnd - tStart) / tSamp;
395  const unsigned n_steps = (unsigned) floor(_n_steps+0.5);
396  assert(fabs(tStart+n_steps*tSamp - tEnd) < 1e-12);
397  const unsigned n_small_steps = (unsigned) floor(tSamp/mDt+0.5);
398  assert(fabs(mDt*n_small_steps - tSamp) < 1e-12);
399 
400  // Initialise solution store
401  OdeSolution solutions;
402  solutions.SetNumberOfTimeSteps(n_steps);
403  solutions.rGetSolutions().push_back(rGetStateVariables());
404  solutions.rGetTimes().push_back(tStart);
405  solutions.SetOdeSystemInformation(this->mpSystemInfo);
406 
407  // Loop over time
408  double curr_time = tStart;
409  for (unsigned i=0; i<n_steps; i++)
410  {
411  for (unsigned j=0; j<n_small_steps; j++)
412  {
413  curr_time = tStart + i*tSamp + j*mDt;
414 
415  // Compute next value of V
416  UpdateTransmembranePotential(curr_time);
417 
418  // Compute other state variables
419  ComputeOneStepExceptVoltage(curr_time);
420 
421  // check gating variables are still in range
423  }
424 
425  // Update solutions
426  solutions.rGetSolutions().push_back(rGetStateVariables());
427  solutions.rGetTimes().push_back(curr_time+mDt);
428  }
429 
430  return solutions;
431  }
432 
442  void ComputeExceptVoltage(double tStart, double tEnd)
443  {
444  // This method iterates over timesteps, calling ComputeExceptVoltage(t) at
445  // each one, to update all state variables except for V, using backward Euler.
446 
447  // Check length of time interval
448  unsigned n_steps = (unsigned)((tEnd - tStart) / mDt + 0.5);
449  assert(fabs(tStart + n_steps*mDt - tEnd) < 1e-12);
450 
451  // Loop over time
452  double curr_time;
453  for (unsigned i=0; i<n_steps; i++)
454  {
455  curr_time = tStart + i*mDt;
456 
457  // Compute other state variables
458  ComputeOneStepExceptVoltage(curr_time);
459 
460 #ifndef NDEBUG
461  // Check gating variables are still in range
463 #endif // NDEBUG
464  }
465  }
466 
474  void SolveAndUpdateState(double tStart, double tEnd)
475  {
476  TimeStepper stepper(tStart, tEnd, mDt);
477 
478  while (!stepper.IsTimeAtEnd())
479  {
480  double time = stepper.GetTime();
481 
482  // Compute next value of V
484 
485  // Compute other state variables
487 
488  // Check gating variables are still in range
490 
491  stepper.AdvanceOneTimeStep();
492  }
493  }
494 
495 private:
496 // LCOV_EXCL_START
504  void EvaluateYDerivatives(double time, const std::vector<double> &rY, std::vector<double> &rDY)
505  {
507  }
508 // LCOV_EXCL_STOP
509 
510 protected:
519  virtual void ComputeOneStepExceptVoltage(double tStart)=0;
520 
528  virtual void UpdateTransmembranePotential(double time)=0;
529 };
530 
531 
532 // Debugging
533 #ifndef NDEBUG
534 #include "OutputFileHandler.hpp"
535 #include <boost/foreach.hpp>
536 template<unsigned SIZE>
537 void DumpJacobianToFile(double time, const double rCurrentGuess[SIZE], double rJacobian[SIZE][SIZE],
538  const std::vector<double>& rY)
539 {
540  OutputFileHandler handler("DumpJacs", false);
541  out_stream p_file = handler.OpenOutputFile("J.txt", std::ios::app);
542  (*p_file) << "At " << time << " " << SIZE << std::endl;
543  (*p_file) << "rY";
544  BOOST_FOREACH(double y, rY)
545  {
546  (*p_file) << " " << y;
547  }
548  (*p_file) << std::endl;
549  (*p_file) << "rCurrentGuess";
550  for (unsigned i=0; i<SIZE; i++)
551  {
552  (*p_file) << " " << rCurrentGuess[i];
553  }
554  (*p_file) << std::endl;
555  (*p_file) << "rJacobian";
556  for (unsigned i=0; i<SIZE; i++)
557  {
558  for (unsigned j=0; j<SIZE; j++)
559  {
560  (*p_file) << " " << rJacobian[i][j];
561  }
562  }
563  (*p_file) << std::endl;
564 }
565 #endif
566 
568 
569 #endif /*ABSTRACTBACKWARDEULERCARDIACCELL_HPP_*/
OdeSolution Compute(double tStart, double tEnd, double tSamp=0.0)
friend class boost::serialization::access
OdeSolution Compute(double tStart, double tEnd, double tSamp=0.0)
#define TEMPLATED_CLASS_IS_ABSTRACT_1_UNSIGNED(T)
std::vector< std::vector< double > > & rGetSolutions()
double GetTime() const
AbstractBackwardEulerCardiacCell(unsigned numberOfStateVariables, unsigned voltageIndex, boost::shared_ptr< AbstractStimulusFunction > pIntracellularStimulus)
void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)
#define NEVER_REACHED
Definition: Exception.hpp:206
void AdvanceOneTimeStep()
bool IsTimeAtEnd() const
void serialize(Archive &archive, const unsigned int version)
void SetOdeSystemInformation(boost::shared_ptr< const AbstractOdeSystemInformation > pOdeSystemInfo)
Definition: OdeSolution.cpp:66
AbstractBackwardEulerCardiacCell(unsigned numberOfStateVariables, unsigned voltageIndex, boost::shared_ptr< AbstractStimulusFunction > pIntracellularStimulus)
virtual void ComputeResidual(double time, const double rCurrentGuess[SIZE], double rResidual[SIZE])=0
virtual void UpdateTransmembranePotential(double time)=0
void SolveAndUpdateState(double tStart, double tEnd)
void ComputeExceptVoltage(double tStart, double tEnd)
std::vector< double > & rGetTimes()
void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)
void serialize(Archive &archive, const unsigned int version)
boost::shared_ptr< AbstractOdeSystemInformation > mpSystemInfo
virtual void ComputeOneStepExceptVoltage(double tStart)=0
void SetNumberOfTimeSteps(unsigned numTimeSteps)
Definition: OdeSolution.cpp:58
virtual void ComputeJacobian(double time, const double rCurrentGuess[SIZE], double rJacobian[SIZE][SIZE])=0