AbstractOneStepIvpOdeSolver.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2010
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 "AbstractOneStepIvpOdeSolver.hpp"
00030 #include "TimeStepper.hpp"
00031 #include <cmath>
00032 
00033 const double smidge = 1e-10;
00034 
00035 OdeSolution AbstractOneStepIvpOdeSolver::Solve(AbstractOdeSystem* pOdeSystem,
00036                                                std::vector<double>& rYValues,
00037                                                double startTime,
00038                                                double endTime,
00039                                                double timeStep,
00040                                                double timeSampling)
00041 {
00042     assert(rYValues.size()==pOdeSystem->GetNumberOfStateVariables());
00043     assert(endTime > startTime);
00044     assert(timeStep > 0.0);
00045     assert(timeSampling >= timeStep);
00046 
00047     mStoppingEventOccurred = false;
00048     if ( pOdeSystem->CalculateStoppingEvent(startTime, rYValues) == true )
00049     {
00050         EXCEPTION("(Solve with sampling) Stopping event is true for initial condition");
00051     }
00052     TimeStepper stepper(startTime, endTime, timeSampling);
00053 
00054     // setup solutions if output is required
00055     OdeSolution solutions;
00056     solutions.SetNumberOfTimeSteps(stepper.EstimateTimeSteps());
00057     solutions.rGetSolutions().push_back(rYValues);
00058     solutions.rGetTimes().push_back(startTime);
00059     solutions.SetOdeSystemInformation(pOdeSystem->GetSystemInformation());
00060 
00061     mWorkingMemory.resize(rYValues.size());
00062 
00063     // Solve the ODE system
00064     while ( !stepper.IsTimeAtEnd() && !mStoppingEventOccurred )
00065     {
00066         InternalSolve(pOdeSystem, rYValues, mWorkingMemory, stepper.GetTime(), stepper.GetNextTime(), timeStep);
00067         stepper.AdvanceOneTimeStep();
00068         // write current solution into solutions
00069         solutions.rGetSolutions().push_back(rYValues);
00070         // Push back new time into the time solution vector
00071         if ( mStoppingEventOccurred )
00072         {
00073             solutions.rGetTimes().push_back(mStoppingTime);
00074         }
00075         else
00076         {
00077             solutions.rGetTimes().push_back(stepper.GetTime());
00078         }
00079     }
00080 
00081     // stepper.EstimateTimeSteps may have been an overestimate...
00082     solutions.SetNumberOfTimeSteps(stepper.GetTotalTimeStepsTaken());
00083     return solutions;
00084 }
00085 
00086 void AbstractOneStepIvpOdeSolver::Solve(AbstractOdeSystem* pOdeSystem,
00087                                         std::vector<double>& rYValues,
00088                                         double startTime,
00089                                         double endTime,
00090                                         double timeStep)
00091 {
00092     assert(rYValues.size()==pOdeSystem->GetNumberOfStateVariables());
00093     assert(endTime > startTime);
00094     assert(timeStep > 0.0);
00095 
00096     mStoppingEventOccurred = false;
00097     if ( pOdeSystem->CalculateStoppingEvent(startTime, rYValues) == true )
00098     {
00099         EXCEPTION("(Solve without sampling) Stopping event is true for initial condition");
00100     }
00101 
00102     // Perhaps resize working memory
00103     mWorkingMemory.resize(rYValues.size());
00104     // And solve...
00105     InternalSolve(pOdeSystem, rYValues, mWorkingMemory, startTime, endTime, timeStep);
00106 }
00107 
00108 void AbstractOneStepIvpOdeSolver::InternalSolve(AbstractOdeSystem* pOdeSystem,
00109                                                 std::vector<double>& rYValues,
00110                                                 std::vector<double>& rWorkingMemory,
00111                                                 double startTime,
00112                                                 double endTime,
00113                                                 double timeStep)
00114 {
00115     TimeStepper stepper(startTime, endTime, timeStep);
00116     // Solve the ODE system
00117 
00118     // Which of our vectors holds the current solution?
00119     // If this is true, it's in rYValues, otherwise it's in rWorkingMemory.
00120     bool curr_is_curr = false;
00121 
00122     // should never get here if this bool has been set to true;
00123     assert(!mStoppingEventOccurred);
00124     while ( !stepper.IsTimeAtEnd() && !mStoppingEventOccurred )
00125     {
00126         curr_is_curr = not curr_is_curr;
00127         // Function that calls the appropriate one-step solver
00128         CalculateNextYValue(pOdeSystem,
00129                             stepper.GetNextTimeStep(),
00130                             stepper.GetTime(),
00131                             curr_is_curr ? rYValues : rWorkingMemory,
00132                             curr_is_curr ? rWorkingMemory : rYValues);
00133         stepper.AdvanceOneTimeStep();
00134         if ( pOdeSystem->CalculateStoppingEvent(stepper.GetTime(),
00135                                                 curr_is_curr ? rWorkingMemory : rYValues) == true )
00136         {
00137             mStoppingTime = stepper.GetTime();
00138             mStoppingEventOccurred = true;
00139         }
00140     }
00141     // Final answer must be in rYValues
00142     if (curr_is_curr)
00143     {
00144         rYValues.assign(rWorkingMemory.begin(), rWorkingMemory.end());
00145     }
00146 }

Generated by  doxygen 1.6.2