CardiacNewtonSolver.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 #ifndef CARDIACNEWTONSOLVER_HPP_
00029 #define CARDIACNEWTONSOLVER_HPP_
00030 
00031 #include <cmath>
00032 #include "IsNan.hpp"
00033 #include "UblasCustomFunctions.hpp"
00034 #include "AbstractBackwardEulerCardiacCell.hpp"
00035 #include "Warnings.hpp"
00036 
00049 template<unsigned SIZE, typename CELLTYPE>
00050 class CardiacNewtonSolver
00051 {
00052 public:
00058     static CardiacNewtonSolver<SIZE, CELLTYPE>* Instance()
00059     {
00060         static CardiacNewtonSolver<SIZE, CELLTYPE> inst;
00061         return &inst;
00062     }
00063 
00071     void Solve(CELLTYPE &rCell,
00072                double time,
00073                double rCurrentGuess[SIZE])
00074     {
00075         unsigned counter = 0;
00076         const double eps = 1e-6; // JonW tolerance
00077 
00078         // check that the initial guess that was given gives a valid residual
00079         rCell.ComputeResidual(time, rCurrentGuess, mResidual.data());
00080         double norm_of_residual = norm_inf(mResidual);
00081         assert(!std::isnan(norm_of_residual));
00082         double norm_of_update = 0.0; //Properly initialised in the loop
00083         do
00084         {
00085             // Calculate Jacobian for current guess
00086             rCell.ComputeJacobian(time, rCurrentGuess, mJacobian);
00087 
00088             // Solve Newton linear system for mUpdate, given mJacobian and mResidual
00089             SolveLinearSystem();
00090 
00091             // Update norm (JonW style)
00092             norm_of_update = norm_inf(mUpdate);
00093 
00094             // Update current guess and recalculate residual
00095             for (unsigned i=0; i<SIZE; i++)
00096             {
00097                 rCurrentGuess[i] -= mUpdate[i];
00098             }
00099             double norm_of_previous_residual = norm_of_residual;
00100             rCell.ComputeResidual(time, rCurrentGuess, mResidual.data());
00101             norm_of_residual = norm_inf(mResidual);
00102             if (norm_of_residual > norm_of_previous_residual && norm_of_update > eps)
00103             {
00104                 //Second part of guard:
00105                 //Note that if norm_of_update < eps (converged) then it's
00106                 //likely that both the residual and the previous residual were
00107                 //close to the root.
00108 
00109                 //Work out where the biggest change in the guess has happened.
00110                 double relative_change_max = 0.0;
00111                 unsigned relative_change_direction = 0;
00112                 for (unsigned i=0; i<SIZE; i++)
00113                 {
00114                     double relative_change = fabs(mUpdate[i]/rCurrentGuess[i]);
00115                     if (relative_change > relative_change_max)
00116                     {
00117                         relative_change_max = relative_change;
00118                         relative_change_direction = i;
00119                     }
00120                 }
00121 
00122                 if (relative_change_max > 1.0)
00123                 {
00124                     //Only walk 0.2 of the way in that direction (put back 0.8)
00125                     rCurrentGuess[relative_change_direction] += 0.8*mUpdate[relative_change_direction];
00126                     rCell.ComputeResidual(time, rCurrentGuess, mResidual.data());
00127                     norm_of_residual = norm_inf(mResidual);
00128                     WARNING("Residual increasing and one direction changing radically - back tracking in that direction");
00129                 }
00130             }
00131             counter++;
00132 
00133             // avoid infinite loops
00134             if (counter > 15)
00135             {
00136 #define COVERAGE_IGNORE
00137                 EXCEPTION("Newton method diverged in CardiacNewtonSolver::Solve()");
00138 #undef COVERAGE_IGNORE
00139             }
00140         }
00141         while (norm_of_update > eps);
00142 
00143 #define COVERAGE_IGNORE
00144 #ifndef NDEBUG
00145         if (norm_of_residual > 2e-10)
00146         { //This line is for correlation - in case we use norm_of_residual as convergence criterion
00147             WARN_ONCE_ONLY("Newton iteration terminated because update vector norm is small, but residual norm is not small.");
00148         }
00149 #endif // NDEBUG
00150 #undef COVERAGE_IGNORE
00151     }
00152 
00153 protected:
00155     CardiacNewtonSolver()
00156     {}
00158     CardiacNewtonSolver(const CardiacNewtonSolver<SIZE, CELLTYPE>&);
00160     CardiacNewtonSolver<SIZE, CELLTYPE>& operator= (const CardiacNewtonSolver<SIZE, CELLTYPE>&);
00161 
00171     void SolveLinearSystem()
00172     {
00173         for (unsigned i=0; i<SIZE; i++)
00174         {
00175             for (unsigned ii=i+1; ii<SIZE; ii++)
00176             {
00177                 double fact = mJacobian[ii][i]/mJacobian[i][i];
00178                 for (unsigned j=i; j<SIZE; j++)
00179                 {
00180                     mJacobian[ii][j] -= fact*mJacobian[i][j];
00181                 }
00182                 mResidual[ii] -= fact*mResidual[i];
00183             }
00184         }
00185         for (unsigned i=SIZE; i-- > 0; )
00186         {
00187             mUpdate[i] = mResidual[i];
00188             for (unsigned j=i+1; j<SIZE; j++)
00189             {
00190                 mUpdate[i] -= mJacobian[i][j]*mUpdate[j];
00191             }
00192             mUpdate[i] /= mJacobian[i][i];
00193         }
00194     }
00195 
00196 private:
00198     c_vector<double, SIZE> mResidual;
00200     double mJacobian[SIZE][SIZE];
00202     c_vector<double, SIZE> mUpdate;
00203 };
00204 
00205 #endif /*CARDIACNEWTONSOLVER_HPP_*/

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