AbstractBackwardEulerCardiacCell.hpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 
00030 #ifndef ABSTRACTBACKWARDEULERCARDIACCELL_HPP_
00031 #define ABSTRACTBACKWARDEULERCARDIACCELL_HPP_
00032 
00033 #include <cassert>
00034 #include <cmath>
00035 
00036 #include "ChasteSerialization.hpp"
00037 #include <boost/serialization/base_object.hpp>
00038 #include "ClassIsAbstract.hpp"
00039 
00040 #include "AbstractCardiacCell.hpp"
00041 #include "Exception.hpp"
00042 #include "PetscTools.hpp"
00043 #include "TimeStepper.hpp"
00044 
00059 template<unsigned SIZE>
00060 class AbstractBackwardEulerCardiacCell : public AbstractCardiacCell
00061 {
00062     private:
00064     friend class boost::serialization::access;
00071     template<class Archive>
00072     void serialize(Archive & archive, const unsigned int version)
00073     {
00074         // This calls serialize on the base class.
00075         archive & boost::serialization::base_object<AbstractCardiacCell>(*this);
00076     }
00077 public:
00078 
00094     AbstractBackwardEulerCardiacCell(
00095         unsigned numberOfStateVariables,
00096         unsigned voltageIndex,
00097         boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus);
00098 
00100     virtual ~AbstractBackwardEulerCardiacCell();
00101 
00109     virtual void ComputeResidual(double time, const double rCurrentGuess[SIZE], double rResidual[SIZE])=0;
00110 
00118     virtual void ComputeJacobian(double time, const double rCurrentGuess[SIZE], double rJacobian[SIZE][SIZE])=0;
00119 
00132     OdeSolution Compute(double tStart, double tEnd, double tSamp=0.0);
00133 
00143     void ComputeExceptVoltage(double tStart, double tEnd);
00144 
00152     void SolveAndUpdateState(double tStart, double tEnd);
00153 
00154 private:
00155 #define COVERAGE_IGNORE
00156 
00163     void EvaluateYDerivatives(double time, const std::vector<double> &rY, std::vector<double> &rDY)
00164     {
00165         NEVER_REACHED;
00166     }
00167 #undef COVERAGE_IGNORE
00168 
00169 protected:
00178     virtual void ComputeOneStepExceptVoltage(double tStart)=0;
00179 
00187     virtual void UpdateTransmembranePotential(double time)=0;
00188 };
00189 
00190 
00191 /*
00192  * NOTE: Explicit instantiation is not used for this class, because the SIZE
00193  * template parameter could take arbitrary values.
00194  */
00195 
00196 
00197 template <unsigned SIZE>
00198 AbstractBackwardEulerCardiacCell<SIZE>::AbstractBackwardEulerCardiacCell(
00199     unsigned numberOfStateVariables,
00200     unsigned voltageIndex,
00201     boost::shared_ptr<AbstractStimulusFunction> pIntracellularStimulus)
00202         : AbstractCardiacCell(boost::shared_ptr<AbstractIvpOdeSolver>(),
00203                               numberOfStateVariables,
00204                               voltageIndex,
00205                               pIntracellularStimulus)
00206 {}
00207 
00208 template <unsigned SIZE>
00209 AbstractBackwardEulerCardiacCell<SIZE>::~AbstractBackwardEulerCardiacCell()
00210 {}
00211 
00212 template <unsigned SIZE>
00213 OdeSolution AbstractBackwardEulerCardiacCell<SIZE>::Compute(double tStart, double tEnd, double tSamp)
00214 {
00215     // In this method, we iterate over timesteps, doing the following for each:
00216     //   - update V using a forward Euler step
00217     //   - call ComputeExceptVoltage(t) to update the remaining state variables
00218     //     using backward Euler
00219 
00220     // Check length of time interval
00221     if (tSamp < mDt)
00222     {
00223         tSamp = mDt;
00224     }
00225     double _n_steps = (tEnd - tStart) / tSamp;
00226     const unsigned n_steps = (unsigned) floor(_n_steps+0.5);
00227     assert(fabs(tStart+n_steps*tSamp - tEnd) < 1e-12);
00228     const unsigned n_small_steps = (unsigned) floor(tSamp/mDt+0.5);
00229     assert(fabs(mDt*n_small_steps - tSamp) < 1e-12);
00230 
00231     // Initialise solution store
00232     OdeSolution solutions;
00233     solutions.SetNumberOfTimeSteps(n_steps);
00234     solutions.rGetSolutions().push_back(rGetStateVariables());
00235     solutions.rGetTimes().push_back(tStart);
00236     solutions.SetOdeSystemInformation(this->mpSystemInfo);
00237 
00238     // Loop over time
00239     double curr_time = tStart;
00240     for (unsigned i=0; i<n_steps; i++)
00241     {
00242         for (unsigned j=0; j<n_small_steps; j++)
00243         {
00244             curr_time = tStart + i*tSamp + j*mDt;
00245 
00246             // Compute next value of V
00247             UpdateTransmembranePotential(curr_time);
00248 
00249             // Compute other state variables
00250             ComputeOneStepExceptVoltage(curr_time);
00251 
00252             // check gating variables are still in range
00253             VerifyStateVariables();
00254         }
00255 
00256         // Update solutions
00257         solutions.rGetSolutions().push_back(rGetStateVariables());
00258         solutions.rGetTimes().push_back(curr_time+mDt);
00259     }
00260 
00261     return solutions;
00262 }
00263 
00264 template <unsigned SIZE>
00265 void AbstractBackwardEulerCardiacCell<SIZE>::ComputeExceptVoltage(double tStart, double tEnd)
00266 {
00267     // This method iterates over timesteps, calling ComputeExceptVoltage(t) at
00268     // each one, to update all state variables except for V, using backward Euler.
00269 
00270     // Check length of time interval
00271     unsigned n_steps = (unsigned)((tEnd - tStart) / mDt + 0.5);
00272     assert(fabs(tStart + n_steps*mDt - tEnd) < 1e-12);
00273 
00274     // Loop over time
00275     double curr_time;
00276     for (unsigned i=0; i<n_steps; i++)
00277     {
00278         curr_time = tStart + i*mDt;
00279 
00280         // Compute other state variables
00281         ComputeOneStepExceptVoltage(curr_time);
00282 
00283 #ifndef NDEBUG
00284         // Check gating variables are still in range
00285         VerifyStateVariables();
00286 #endif // NDEBUG
00287     }
00288 }
00289 
00290 template<unsigned SIZE>
00291 void AbstractBackwardEulerCardiacCell<SIZE>::SolveAndUpdateState(double tStart, double tEnd)
00292 {
00293     TimeStepper stepper(tStart, tEnd, mDt);
00294 
00295     while(!stepper.IsTimeAtEnd())
00296     {
00297         double time = stepper.GetTime();
00298 
00299         // Compute next value of V
00300         UpdateTransmembranePotential(time);
00301 
00302         // Compute other state variables
00303         ComputeOneStepExceptVoltage(time);
00304 
00305         // check gating variables are still in range
00306         VerifyStateVariables();
00307 
00308         stepper.AdvanceOneTimeStep();
00309     }
00310 }
00311 
00312 
00313 TEMPLATED_CLASS_IS_ABSTRACT_1_UNSIGNED(AbstractBackwardEulerCardiacCell)
00314 
00315 #endif /*ABSTRACTBACKWARDEULERCARDIACCELL_HPP_*/

Generated on Mon Apr 18 11:35:28 2011 for Chaste by  doxygen 1.5.5