NhsSystemWithImplicitSolver.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2009
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 #include "NhsSystemWithImplicitSolver.hpp"
00029 #include "TimeStepper.hpp"
00030 #include <cmath>
00031 #include "LogFile.hpp"
00032 #include <iostream>
00033 /*
00034  *=========================== PRIVATE METHODS ==============================
00035  */
00036 
00037 void NhsSystemWithImplicitSolver::ImplicitSolveForActiveTension()
00038 {
00039     // solve a 1d nonlinear problem for the active tension
00040 
00041     // current active tension
00042     double current_active_tension = mActiveTensionInitialGuess;
00043 
00044     // see what the residual is
00045     double residual = CalcActiveTensionResidual(current_active_tension);
00046 
00047     // solve using Newton's method, no damping. Stop if num iterations
00048     // reaches 15 (very conservative)
00049     unsigned counter = 0;
00050     std::vector<double> old_residuals(15);
00051     while ((fabs(residual)>mTolerance) && (counter++<15))
00052     {
00053         // numerically approximate the jacobian
00054         double h = std::max(fabs(current_active_tension/100),1e-8);
00055 
00056         double jac = (CalcActiveTensionResidual(current_active_tension+h) - CalcActiveTensionResidual(current_active_tension))/h;
00057 
00058         current_active_tension -= residual/jac;
00059         residual = CalcActiveTensionResidual(current_active_tension);
00060         old_residuals.push_back(residual);
00061     }
00062     //assert(counter<15);
00063     if(counter >= 15) // not sure what to do here, after having got a case where resid stagnated on +/- 5.500502e-10
00064     {
00065         #define COVERAGE_IGNORE
00066         LOG(1, "\nWARNING in NhsSystemWithImplicitSolver::ImplicitSolveForActiveTension(), counter="<<counter<<",resids=");
00067         for (unsigned i=0; i<old_residuals.size(); i++)
00068         {
00069             LOG(1,old_residuals[i]);
00070         }
00071         LOG(1,"Final residual = " << residual);
00072         if(residual > 100*mTolerance)
00073         {
00074             LOG(1,"Residual > 100*mTol, throwing exception\n");
00075             EXCEPTION("NhsSystemWithImplicitSolver::ImplicitSolveForActiveTension() failed to converge");
00076         }
00077         #undef COVERAGE_IGNORE
00078     }
00079 
00080 
00081     // save the active tension initial guess for next time round
00082     mActiveTensionInitialGuess = current_active_tension;
00083 
00087     mActiveTensionSolution = current_active_tension;
00088 }
00089 
00090 double NhsSystemWithImplicitSolver::CalcActiveTensionResidual(double activeTensionGuess)
00091 {
00092     // to calculate the active tension residual, we solve use the current active tension,
00093     // solve for Ca_trop implicitly, then z implicitly, then Q implicitly, then see
00094     // what the resulting active tension was
00095 
00096     // solve for new Ca_trop implicitly
00097     double new_Ca_Trop = ImplicitSolveForCaTrop(activeTensionGuess);
00098 
00099     // store the current Ca_trop solution, in case it turns out to be the current one
00100     mTempStoredStateVariables[0] = new_Ca_Trop;
00101 
00102     // solve for new z implicitly (or with a mixed implicit-explicit scheme, if asked for)
00103     double new_z;
00104     if(!mUseImplicitExplicitSolveForZ)
00105     {
00106         new_z = ImplicitSolveForZ(new_Ca_Trop);
00107     }
00108     else
00109     {
00110         new_z = ImplicitExplicitSolveForZ(new_Ca_Trop);
00111     }
00112 
00113     // store the current z solution, in case it turns out to be the current one
00114     mTempStoredStateVariables[1] = new_z;
00115 
00116     // get the new T0
00117     double new_T0 = CalculateT0(new_z);
00118 
00119     // solve for the new Qi's (and therefore Q) implicitly (note Q1, Q2, Q3 are stored in
00120     // this method
00121     double new_Q = ImplicitSolveForQ();                        // <-- SHOULD BE MOVED UP!!
00122     //double new_Q = QuasiStaticSolveForQ();
00123 
00124 
00125     // compute the new active tension and return the difference
00126     double new_active_tension = 0;
00127     if(new_Q > 0)
00128     {
00129         new_active_tension = new_T0*(1+(2+mA)*new_Q)/(1+new_Q);
00130     }
00131     else
00132     {
00133         #define COVERAGE_IGNORE // not quite sure how to cover this, of if it is worth it
00134         new_active_tension = new_T0*(1+mA*new_Q)/(1-new_Q);
00135         #undef COVERAGE_IGNORE
00136     }
00137 
00138     return new_active_tension - activeTensionGuess;
00139 }
00140 
00141 
00142 // Assume the active tension is known and solve for the Ca_trop at the next time
00143 // implicitly using backward euler. This can be done directly as the rhs is linear
00144 // in Ca_trop
00145 double NhsSystemWithImplicitSolver::ImplicitSolveForCaTrop(double newActiveTension)
00146 {
00147     double old_Ca_trop = mCurrentStateVars[0];
00148 
00149     double numer = old_Ca_trop + mDt * mKon * mCalciumI  * mCalciumTroponinMax;
00150     double denom = 1 + mDt*mKon*mCalciumI + mDt*mKrefoff*(1-newActiveTension/(mGamma*mTref));
00151 
00152     return numer/denom;
00153 }
00154 
00155 // Assume the Ca_trop is known and solve for the z at the next time
00156 // implicitly using backward euler. Uses Newton's method
00157 double NhsSystemWithImplicitSolver::ImplicitSolveForZ(double newCaTrop)
00158 {
00159     double current_z = mCurrentStateVars[1];
00160     double residual = CalcZResidual(current_z, newCaTrop);
00161     unsigned counter = 0;
00162 
00163     while ((fabs(residual)>mTolerance) && (counter++<15))
00164     {
00165         double h = std::max(fabs(current_z/100),1e-8);
00166         double jac = (CalcZResidual(current_z + h, newCaTrop) - CalcZResidual(current_z, newCaTrop));
00167         jac/= h;
00168 
00169         current_z -= residual/jac;
00170         residual = CalcZResidual(current_z, newCaTrop);
00171     }
00172     assert(counter<15);
00173 
00174     return current_z;
00175 }
00176 
00177 
00178 double NhsSystemWithImplicitSolver::CalcZResidual(double z, double newCaTrop)
00179 {
00180     double ca_trop_to_ca_trop50_ratio_to_n = pow(newCaTrop/mCalciumTrop50, mN);
00181     double dzdt =  mAlpha0 * ca_trop_to_ca_trop50_ratio_to_n * (1-z)
00182                  - mAlphaR1 * z
00183                  - mAlphaR2 * pow(z,mNr) / (pow(z,mNr) + pow(mKZ,mNr));
00184 
00185     return z - mCurrentStateVars[1] - mDt*dzdt;
00186 }
00187 
00188 
00189 // Solve for z semi-implicitly instead of fully implicitly. If we assume we know
00190 // Ca_trop solving for z is a 1d nonlinear problem. Call this to treat the problem
00191 // implicitly in the linear terms on the rhs of dzdt (the (1-z) and (z) terms), and
00192 // explicitly in the nonlinear term (the z^nr/(z^nr + K^nr) term. This means the
00193 // problem can be solved directly and no Newton iterations are needed.
00194 double NhsSystemWithImplicitSolver::ImplicitExplicitSolveForZ(double newCaTrop)
00195 {
00196     double ca_trop_to_ca_trop50_ratio_to_n = pow(newCaTrop/mCalciumTrop50, mN);
00197     double old_z = mCurrentStateVars[1];
00198 
00199     double numer =   old_z + mDt * mAlpha0 * ca_trop_to_ca_trop50_ratio_to_n
00200                    - mDt * mAlphaR2 * pow(old_z,mNr) / (pow(old_z,mNr) + pow(mKZ,mNr));
00201     double denom =   1 + mDt * mAlpha0 * ca_trop_to_ca_trop50_ratio_to_n + mAlphaR1 *mDt;
00202 
00203     return numer/denom;
00204 }
00205 
00206 
00207 // Solve for Q1, Q2, Q3 implicitly and directly
00208 double NhsSystemWithImplicitSolver::ImplicitSolveForQ()
00209 {
00210     double old_Q1 = mCurrentStateVars[2];
00211     double old_Q2 = mCurrentStateVars[3];
00212     double old_Q3 = mCurrentStateVars[4];
00213 
00214     double new_Q1 = (old_Q1 + mDt*mA1*mDLambdaDt)/(1 + mAlpha1*mDt);
00215     double new_Q2 = (old_Q2 + mDt*mA2*mDLambdaDt)/(1 + mAlpha2*mDt);
00216     double new_Q3 = (old_Q3 + mDt*mA3*mDLambdaDt)/(1 + mAlpha3*mDt);
00217 
00218     mTempStoredStateVariables[2] = new_Q1;
00219     mTempStoredStateVariables[3] = new_Q2;
00220     mTempStoredStateVariables[4] = new_Q3;
00221 
00222     return new_Q1 + new_Q2 + new_Q3;
00223 }
00224 
00225 
00226 //double NhsSystemWithImplicitSolver::QuasiStaticSolveForQ()
00227 //{
00228 //    double new_Q1 = mA1*mDLambdaDt/mAlpha1;
00229 //    double new_Q2 = mA2*mDLambdaDt/mAlpha2;
00230 //    double new_Q3 = mA3*mDLambdaDt/mAlpha3;
00231 //
00232 //    mTempStoredStateVariables[2] = new_Q1;
00233 //    mTempStoredStateVariables[3] = new_Q2;
00234 //    mTempStoredStateVariables[4] = new_Q3;
00235 //
00236 //    return new_Q1 + new_Q2 + new_Q3;
00237 //}
00238 
00239 
00240 
00241 /*
00242  *=========================== PUBLIC METHODS ==============================
00243  */
00244 
00245 NhsSystemWithImplicitSolver::NhsSystemWithImplicitSolver()
00246     : NhsCellularMechanicsOdeSystem()
00247 {
00248     mTempStoredStateVariables.resize(GetNumberOfStateVariables());
00249     mActiveTensionInitialGuess = GetActiveTension();
00250     mUseImplicitExplicitSolveForZ = false;
00251 
00252     mActiveTensionSolution = NhsCellularMechanicsOdeSystem::GetActiveTension();
00253 }
00254 
00255 void NhsSystemWithImplicitSolver::SetActiveTensionInitialGuess(double activeTensionInitialGuess)
00256 {
00257     mActiveTensionInitialGuess = activeTensionInitialGuess;
00258 }
00259 
00260 
00261 void NhsSystemWithImplicitSolver::SolveDoNotUpdate(double startTime,
00262                                                    double endTime,
00263                                                    double timestep)
00264 {
00265     assert(startTime < endTime);
00266 
00267     mDt = timestep;
00268 
00269     // the implicit routines look in mCurrentStateVars for the current values
00270     mCurrentStateVars = mStateVariables;
00271 
00272     // loop in time
00273     TimeStepper stepper(startTime, endTime, timestep);
00274     while ( !stepper.IsTimeAtEnd() )
00275     {
00276         ImplicitSolveForActiveTension();
00277         mCurrentStateVars = mTempStoredStateVariables;
00278 
00279         stepper.AdvanceOneTimeStep();
00280     }
00281 
00282     // note: state variables NOT updated.
00283 }
00284 
00285 void NhsSystemWithImplicitSolver::UpdateStateVariables()
00286 {
00287     for(unsigned i=0; i<mStateVariables.size(); i++)
00288     {
00289         mStateVariables[i] = mCurrentStateVars[i];
00290     }
00291 }
00292 
00293 void NhsSystemWithImplicitSolver::UseImplicitExplicitSolveForZ(bool useImplicitExplicitSolveForZ)
00294 {
00295     mUseImplicitExplicitSolveForZ = useImplicitExplicitSolveForZ;
00296 }
00297 
00298 double NhsSystemWithImplicitSolver::GetActiveTensionAtNextTime()
00299 {
00300     return mActiveTensionSolution;
00301 }
00302 

Generated on Wed Mar 18 12:51:51 2009 for Chaste by  doxygen 1.5.5