NhsModelWithBackwardSolver.cpp

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 #include "NhsModelWithBackwardSolver.hpp"
00030 #include <iostream>
00031 #include <cmath>
00032 #include "LogFile.hpp"
00033 #include "Exception.hpp"
00034 #include "TimeStepper.hpp"
00035 
00036 const double NhsModelWithBackwardSolver::mTolerance = 1e-10;
00037 
00038 double NhsModelWithBackwardSolver::ImplicitSolveForQ()
00039 {
00040     mTemporaryStateVariables[2] = (mTemporaryStateVariables[2] + mDt*mA1*mDLambdaDt)/(1 + mAlpha1*mDt);
00041     mTemporaryStateVariables[3] = (mTemporaryStateVariables[3] + mDt*mA2*mDLambdaDt)/(1 + mAlpha2*mDt);
00042     mTemporaryStateVariables[4] = (mTemporaryStateVariables[4] + mDt*mA3*mDLambdaDt)/(1 + mAlpha3*mDt);
00043 
00044     return mTemporaryStateVariables[2] + mTemporaryStateVariables[3] + mTemporaryStateVariables[4];
00045 }
00046 
00047 void NhsModelWithBackwardSolver::CalculateCaTropAndZDerivatives(double calciumTroponin, double z, double Q,
00048                                                                 double& dCaTrop, double& dz)
00049 {
00050 //As in straight Nhs, we don't cover the exception code
00051 #define COVERAGE_IGNORE
00052     if(calciumTroponin < 0)
00053     {
00054         EXCEPTION("CalciumTrop concentration went negative");
00055     }
00056     if(z<0)
00057     {
00058         EXCEPTION("z went negative");
00059     }
00060     if(z>1)
00061     {
00062         EXCEPTION("z became greater than 1");
00063     }
00064 #undef COVERAGE_IGNORE
00065 
00066     double T0 = CalculateT0(z);
00067 
00068     double Ta;
00069     if(Q>0)
00070     {
00071         Ta = T0*(1+(2+mA)*Q)/(1+Q);
00072     }
00073     else
00074     {
00075         Ta = T0*(1+mA*Q)/(1-Q);
00076     }
00077 
00078     dCaTrop =   mKon * mCalciumI * ( mCalciumTroponinMax - calciumTroponin)
00079              - mKrefoff * (1-Ta/(mGamma*mTref)) * calciumTroponin;
00080 
00081     double ca_trop_to_ca_trop50_ratio_to_n = pow(calciumTroponin/mCalciumTrop50, mN);
00082 
00083     dz =   mAlpha0 * ca_trop_to_ca_trop50_ratio_to_n * (1-z)
00084          - mAlphaR1 * z
00085          - mAlphaR2 * pow(z,mNr) / (pow(z,mNr) + pow(mKZ,mNr));
00086 }
00087 
00088 
00089 
00090 void NhsModelWithBackwardSolver::CalculateBackwardEulerResidual(double calciumTroponin, double z, double Q,
00091                                                                 double& residualComponent1, double& residualComponent2)
00092 {
00093     double dcatrop;
00094     double dz;
00095     CalculateCaTropAndZDerivatives(calciumTroponin,z,Q,dcatrop,dz);
00096 
00097     residualComponent1 = calciumTroponin - mDt*dcatrop - mTemporaryStateVariables[0];
00098     residualComponent2 = z - mDt*dz - mTemporaryStateVariables[1];
00099 }
00100 
00101 
00102 
00103 NhsModelWithBackwardSolver::NhsModelWithBackwardSolver()
00104 {
00105     mTemporaryStateVariables.resize(5);
00106 }
00107 
00108 
00109 
00110 void NhsModelWithBackwardSolver::RunDoNotUpdate(double startTime, double endTime, double timestep)
00111 {
00112     assert(startTime < endTime);
00113 
00114     mDt = timestep;
00115 
00116     mTemporaryStateVariables = mStateVariables;
00117 
00118     // loop in time
00119     TimeStepper stepper(startTime, endTime, timestep);
00120 
00121     while ( !stepper.IsTimeAtEnd() )
00122     {
00124         // Q1,Q2,Q3 using backward euler can solved straightaway
00126         double new_Q = ImplicitSolveForQ();
00127 
00129         // Solve the 2D nonlinear problem for Backward Euler Ca_trop and z
00131 
00132         // see what the residual is
00133         double catrop_guess = mTemporaryStateVariables[0];
00134         double z_guess = mTemporaryStateVariables[1];
00135         double f1,f2; // f=[f1,f2]=residual
00136 
00137         CalculateBackwardEulerResidual(catrop_guess, z_guess, new_Q, f1, f2);
00138         double norm_resid = sqrt(f1*f1+f2*f2);
00139 
00140         // solve using Newton's method, no damping. Stop if num iterations
00141         // reaches 15 (very conservative)
00142         unsigned counter = 0;
00143         while ((norm_resid>mTolerance) && (counter++<15))
00144         {
00145             // numerically approximate the jacobian J
00146             double j11,j12,j21,j22; // J = [j11, j12; j21 j22]
00147             double temp1,temp2;
00148 
00149             double h = std::max(fabs(catrop_guess/100),1e-8);
00150             CalculateBackwardEulerResidual(catrop_guess+h, z_guess, new_Q, temp1, temp2);
00151             j11 = (temp1-f1)/h;
00152             j21 = (temp2-f2)/h;
00153 
00154             h = std::max(fabs(z_guess/100),1e-8);
00155             CalculateBackwardEulerResidual(catrop_guess, z_guess+h, new_Q, temp1, temp2);
00156             j12 = (temp1-f1)/h;
00157             j22 = (temp2-f2)/h;
00158 
00159             // compute u = J^{-1} f (exactly, as a 2D problem)
00160             double one_over_det = 1.0/(j11*j22-j12*j21);
00161             double u1 = one_over_det*(j22*f1  - j12*f2);
00162             double u2 = one_over_det*(-j21*f1 + j11*f2);
00163 
00164             catrop_guess -= u1;
00165             z_guess -= u2;
00166 
00167             CalculateBackwardEulerResidual(catrop_guess, z_guess, new_Q, f1, f2);
00168             norm_resid = sqrt(f1*f1+f2*f2);
00169         }
00170         assert(counter<15); // if this fails, see corresponding code in old NhsModelWithImplicitSolver
00171 
00172         mTemporaryStateVariables[0] = catrop_guess;
00173         mTemporaryStateVariables[1] = z_guess;
00174 
00175         stepper.AdvanceOneTimeStep();
00176     }
00177 }
00178 
00179 double NhsModelWithBackwardSolver::GetNextActiveTension()
00180 {
00181     double T0 = CalculateT0(mTemporaryStateVariables[1]);
00182     double Q = mTemporaryStateVariables[2]+mTemporaryStateVariables[3]+mTemporaryStateVariables[4];
00183 
00184     if(Q>0)
00185     {
00186         return T0*(1+(2+mA)*Q)/(1+Q);
00187     }
00188     else
00189     {
00190         return T0*(1+mA*Q)/(1-Q);
00191     }
00192 }
00193 
00194 void NhsModelWithBackwardSolver::RunAndUpdate(double startTime, double endTime, double timestep)
00195 {
00196     RunDoNotUpdate(startTime, endTime, timestep);
00197     UpdateStateVariables();
00198 }
00199 

Generated on Tue May 31 14:31:41 2011 for Chaste by  doxygen 1.5.5