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

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