Chaste  Release::2018.1
NhsModelWithBackwardSolver.cpp
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 #include "NhsModelWithBackwardSolver.hpp"
37 #include <iostream>
38 #include <cmath>
39 #include "LogFile.hpp"
40 #include "Exception.hpp"
41 #include "TimeStepper.hpp"
42 
43 const double NhsModelWithBackwardSolver::mTolerance = 1e-10;
44 
46 {
50 
52 }
53 
54 void NhsModelWithBackwardSolver::CalculateCaTropAndZDerivatives(double calciumTroponin, double z, double Q,
55  double& dCaTrop, double& dz)
56 {
57 //As in straight Nhs, we don't cover the exception code
58 // LCOV_EXCL_START
59  if (calciumTroponin < 0)
60  {
61  EXCEPTION("CalciumTrop concentration went negative");
62  }
63  if (z<0)
64  {
65  EXCEPTION("z went negative");
66  }
67  if (z>1)
68  {
69  EXCEPTION("z became greater than 1");
70  }
71 // LCOV_EXCL_STOP
72 
73  double T0 = CalculateT0(z);
74 
75  double Ta;
76  if (Q>0)
77  {
78  Ta = T0*(1+(2+mA)*Q)/(1+Q);
79  }
80  else
81  {
82  Ta = T0*(1+mA*Q)/(1-Q);
83  }
84 
85  dCaTrop = mKon * mCalciumI * ( mCalciumTroponinMax - calciumTroponin)
86  - mKrefoff * (1-Ta/(mGamma*mTref)) * calciumTroponin;
87 
88  double ca_trop_to_ca_trop50_ratio_to_n = pow(calciumTroponin/mCalciumTrop50, mN);
89 
90  dz = mAlpha0 * ca_trop_to_ca_trop50_ratio_to_n * (1-z)
91  - mAlphaR1 * z
92  - mAlphaR2 * pow(z,mNr) / (pow(z,mNr) + pow(mKZ,mNr));
93 }
94 
95 void NhsModelWithBackwardSolver::CalculateBackwardEulerResidual(double calciumTroponin, double z, double Q,
96  double& residualComponent1, double& residualComponent2)
97 {
98  double dcatrop;
99  double dz;
100  CalculateCaTropAndZDerivatives(calciumTroponin,z,Q,dcatrop,dz);
101 
102  residualComponent1 = calciumTroponin - mDt*dcatrop - mTemporaryStateVariables[0];
103  residualComponent2 = z - mDt*dz - mTemporaryStateVariables[1];
104 }
105 
107 {
108  mTemporaryStateVariables.resize(5);
109 }
110 
111 void NhsModelWithBackwardSolver::RunDoNotUpdate(double startTime, double endTime, double timestep)
112 {
113  assert(startTime < endTime);
114 
115  mDt = timestep;
116 
118 
119  // loop in time
120  TimeStepper stepper(startTime, endTime, timestep);
121 
122  while ( !stepper.IsTimeAtEnd() )
123  {
125  // Q1,Q2,Q3 using backward euler can solved straightaway
127  double new_Q = ImplicitSolveForQ();
128 
130  // Solve the 2D nonlinear problem for Backward Euler Ca_trop and z
132 
133  // see what the residual is
134  double catrop_guess = mTemporaryStateVariables[0];
135  double z_guess = mTemporaryStateVariables[1];
136  double f1,f2; // f=[f1,f2]=residual
137 
138  CalculateBackwardEulerResidual(catrop_guess, z_guess, new_Q, f1, f2);
139  double norm_resid = sqrt(f1*f1+f2*f2);
140 
141  // solve using Newton's method, no damping. Stop if num iterations
142  // reaches 15 (very conservative)
143  unsigned counter = 0;
144  while ((norm_resid>mTolerance) && (counter++<15))
145  {
146  // numerically approximate the jacobian J
147  double j11,j12,j21,j22; // J = [j11, j12; j21 j22]
148  double temp1,temp2;
149 
150  double h = std::max(fabs(catrop_guess/100),1e-8);
151  CalculateBackwardEulerResidual(catrop_guess+h, z_guess, new_Q, temp1, temp2);
152  j11 = (temp1-f1)/h;
153  j21 = (temp2-f2)/h;
154 
155  h = std::max(fabs(z_guess/100),1e-8);
156  CalculateBackwardEulerResidual(catrop_guess, z_guess+h, new_Q, temp1, temp2);
157  j12 = (temp1-f1)/h;
158  j22 = (temp2-f2)/h;
159 
160  // compute u = J^{-1} f (exactly, as a 2D problem)
161  double one_over_det = 1.0/(j11*j22-j12*j21);
162  double u1 = one_over_det*(j22*f1 - j12*f2);
163  double u2 = one_over_det*(-j21*f1 + j11*f2);
164 
165  catrop_guess -= u1;
166  z_guess -= u2;
167 
168  CalculateBackwardEulerResidual(catrop_guess, z_guess, new_Q, f1, f2);
169  norm_resid = sqrt(f1*f1+f2*f2);
170  }
171  assert(counter<15); // if this fails, see corresponding code in old NhsModelWithImplicitSolver
172 
173  mTemporaryStateVariables[0] = catrop_guess;
174  mTemporaryStateVariables[1] = z_guess;
175 
176  stepper.AdvanceOneTimeStep();
177  }
178 }
179 
181 {
182  double T0 = CalculateT0(mTemporaryStateVariables[1]);
184 
185  if (Q>0)
186  {
187  return T0*(1+(2+mA)*Q)/(1+Q);
188  }
189  else
190  {
191  return T0*(1+mA*Q)/(1-Q);
192  }
193 }
194 
195 void NhsModelWithBackwardSolver::RunAndUpdate(double startTime, double endTime, double timestep)
196 {
197  RunDoNotUpdate(startTime, endTime, timestep);
199 }
200 
void CalculateBackwardEulerResidual(double calciumTroponin, double z, double Q, double &residualComponent1, double &residualComponent2)
void CalculateCaTropAndZDerivatives(double calciumTroponin, double z, double Q, double &dCaTrop, double &dz)
double CalculateT0(double z)
static const double mTref
#define EXCEPTION(message)
Definition: Exception.hpp:143
static const double mKZ
static const unsigned mNr
void AdvanceOneTimeStep()
static const double mAlphaR1
bool IsTimeAtEnd() const
static const double mAlpha0
static const unsigned mN
static const double mAlpha3
static const double mCalciumTroponinMax
void RunAndUpdate(double startTime, double endTime, double timestep)
static const double mA
void RunDoNotUpdate(double startTime, double endTime, double timestep)
static const double mA2
static const double mA3
static const double mAlpha2
static const double mGamma
static const double mKon
static const double mAlpha1
static const double mAlphaR2
static const double mKrefoff
static const double mA1