Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
AbstractBackwardEulerCardiacCell.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
37#ifndef ABSTRACTBACKWARDEULERCARDIACCELL_HPP_
38#define ABSTRACTBACKWARDEULERCARDIACCELL_HPP_
39
40#include <cassert>
41#include <cmath>
42
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
66template<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 }
84public:
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
161private:
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
176protected:
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
207template <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
218template <unsigned SIZE>
221
222template <unsigned SIZE>
223OdeSolution 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
274template <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
300template<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
329template<>
331{
332private:
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
348public:
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
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
417
418 // Compute other state variables
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
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
495private:
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
510protected:
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>
536template<unsigned SIZE>
537void 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_*/
#define TEMPLATED_CLASS_IS_ABSTRACT_1_UNSIGNED(T)
#define NEVER_REACHED
AbstractBackwardEulerCardiacCell(unsigned numberOfStateVariables, unsigned voltageIndex, boost::shared_ptr< AbstractStimulusFunction > pIntracellularStimulus)
OdeSolution Compute(double tStart, double tEnd, double tSamp=0.0)
void serialize(Archive &archive, const unsigned int version)
void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)
virtual void ComputeOneStepExceptVoltage(double tStart)=0
virtual void UpdateTransmembranePotential(double time)=0
void ComputeExceptVoltage(double tStart, double tEnd)
virtual void ComputeJacobian(double time, const double rCurrentGuess[SIZE], double rJacobian[SIZE][SIZE])=0
void EvaluateYDerivatives(double time, const std::vector< double > &rY, std::vector< double > &rDY)
virtual void ComputeOneStepExceptVoltage(double tStart)=0
void serialize(Archive &archive, const unsigned int version)
virtual void ComputeResidual(double time, const double rCurrentGuess[SIZE], double rResidual[SIZE])=0
void SolveAndUpdateState(double tStart, double tEnd)
OdeSolution Compute(double tStart, double tEnd, double tSamp=0.0)
virtual void UpdateTransmembranePotential(double time)=0
AbstractBackwardEulerCardiacCell(unsigned numberOfStateVariables, unsigned voltageIndex, boost::shared_ptr< AbstractStimulusFunction > pIntracellularStimulus)
boost::shared_ptr< AbstractOdeSystemInformation > mpSystemInfo
void SetNumberOfTimeSteps(unsigned numTimeSteps)
std::vector< std::vector< double > > & rGetSolutions()
std::vector< double > & rGetTimes()
void SetOdeSystemInformation(boost::shared_ptr< const AbstractOdeSystemInformation > pOdeSystemInfo)
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
bool IsTimeAtEnd() const
double GetTime() const
void AdvanceOneTimeStep()