CvodeAdaptor.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 #ifdef CHASTE_CVODE
00030 
00031 #include "CvodeAdaptor.hpp"
00032 
00033 #include "Exception.hpp"
00034 #include "TimeStepper.hpp"
00035 #include "VectorHelperFunctions.hpp"
00036 
00037 #include <iostream>
00038 #include <sstream>
00039 
00040 // CVODE headers
00041 #include <sundials/sundials_nvector.h>
00042 #include <cvode/cvode_dense.h>
00043 
00044 
00061 int CvodeRhsAdaptor(realtype t, N_Vector y, N_Vector ydot, void* pData)
00062 {
00063     assert(pData != NULL);
00064     CvodeData* p_data = (CvodeData*) pData;
00065     // Get y, ydot into std::vector<>s
00066     static std::vector<realtype> ydot_vec;
00067     CopyToStdVector(y, *p_data->pY);
00068     CopyToStdVector(ydot, ydot_vec);
00069     // Call our function
00070     try
00071     {
00072         p_data->pSystem->EvaluateYDerivatives(t, *(p_data->pY), ydot_vec);
00073     }
00074     catch (const Exception &e)
00075     {
00076         std::cerr << "CVODE RHS Exception: " << e.GetMessage() << std::endl << std::flush;
00077         return -1;
00078     }
00079     // Copy derivative back
00080     CopyFromStdVector(ydot_vec, ydot);
00081     return 0;
00082 }
00083 
00104 int CvodeRootAdaptor(realtype t, N_Vector y, realtype* pGOut, void* pData)
00105 {
00106     assert(pData != NULL);
00107     CvodeData* p_data = (CvodeData*) pData;
00108     // Get y into a std::vector
00109     CopyToStdVector(y, *p_data->pY);
00110     // Call our function
00111     try
00112     {
00113         *pGOut = p_data->pSystem->CalculateRootFunction(t, *p_data->pY);
00114     }
00115     catch (const Exception &e)
00116     {
00117         std::cerr << "CVODE Root Exception: " << e.GetMessage() << std::endl << std::flush;
00118         return -1;
00119     }
00120     return 0;
00121 }
00122 
00123 // /**
00124 //  * Jacobian computation adaptor function.
00125 //  *
00126 //  * If solving an AbstractOdeSystemWithAnalyticJacobian, this function
00127 //  * can be used to allow CVODE to compute the Jacobian analytically.
00128 //  *
00129 //  * Note to self: can test using pSystem->GetUseAnalyticJacobian().
00130 //  */
00131 // int CvodeDenseJacobianAdaptor(long int numberOfStateVariables, DenseMat J,
00132 //                               realtype t, N_Vector y, N_Vector fy,
00133 //                               void* pData,
00134 //                               N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
00135 // {
00136 //     AbstractOdeSystemWithAnalyticJacobian* pSystem
00137 //         = (AbstractOdeSystemWithAnalyticJacobian*) pData;
00138 //     // Get the current time step
00139 //     double dt;
00140 //     CVodeGetCurrentStep(CvodeMem, &dt);
00141 //     // Get std::vector<> for y and double** for J
00142 //     std::vector<realtype>& y_vec = *NV_DATA_STL(y);
00143 //     double** ppJ = J->data; // organised column-wise: J_{i,j} = ppJ[j][i]
00144 //     // Call our function
00145 //     try
00146 //     {
00147 //         pSystem->AnalyticJacobian(y_vec, ppJ, t, dt);
00148 //     }
00149 //     catch (const Exception &e)
00150 //     {
00151 //         std::cerr << "CVODE J Exception: " << e.GetMessage() << std::endl << std::flush;
00152 //         return -1;
00153 //     }
00154 //     // Update J (if needed)
00155 //     return 0;
00156 // }
00157 
00158 
00159 
00160 void CvodeErrorHandler(int errorCode, const char *module, const char *function,
00161                        char *message, void* pData)
00162 {
00163     std::stringstream err;
00164     err << "CVODE Error " << errorCode << " in module " << module
00165         << " function " << function << ": " << message;
00166     std::cerr << "*" << err.str() << std::endl << std::flush;
00167     // Throwing the exception here causes termination on Maverick (g++ 4.4)
00168     //EXCEPTION(err.str());
00169 }
00170 
00171 
00172 
00173 void CvodeAdaptor::SetupCvode(AbstractOdeSystem* pOdeSystem,
00174                               std::vector<double>& rInitialY,
00175                               double startTime, double maxStep)
00176 {
00177     assert(rInitialY.size() == pOdeSystem->GetNumberOfStateVariables());
00178     assert(maxStep > 0.0);
00179 
00180     mInitialValues = N_VMake_Serial(rInitialY.size(), &(rInitialY[0]));
00181     assert(NV_DATA_S(mInitialValues) == &(rInitialY[0]));
00182     assert(!NV_OWN_DATA_S(mInitialValues));
00183 //    std::cout << " Initial values: "; N_VPrint_Serial(mInitialValues);
00184 //    std::cout << " Rtol: " << mRelTol << ", Atol: " << mAbsTol << std::endl;
00185 //    std::cout << " Start: " << startTime << " max dt=" << maxStep << std::endl << std::flush;
00186 
00187     mpCvodeMem = CVodeCreate(CV_BDF, CV_NEWTON);
00188     if (mpCvodeMem == NULL) EXCEPTION("Failed to SetupCvode CVODE");
00189     // Set error handler
00190     CVodeSetErrHandlerFn(mpCvodeMem, CvodeErrorHandler, NULL);
00191     // Set the user data
00192     mData.pSystem = pOdeSystem;
00193     mData.pY = &rInitialY;
00194 #if CHASTE_SUNDIALS_VERSION >= 20400
00195     CVodeSetUserData(mpCvodeMem, (void*)(&mData));
00196 #else
00197     CVodeSetFdata(mpCvodeMem, (void*)(&mData));
00198 #endif
00199     // Setup CVODE
00200 #if CHASTE_SUNDIALS_VERSION >= 20400
00201     CVodeInit(mpCvodeMem, CvodeRhsAdaptor, startTime, mInitialValues); 
00202     CVodeSStolerances(mpCvodeMem, mRelTol, mAbsTol);
00203 #else
00204     CVodeMalloc(mpCvodeMem, CvodeRhsAdaptor, startTime, mInitialValues,
00205                 CV_SS, mRelTol, &mAbsTol);
00206 #endif
00207     CVodeSetMaxStep(mpCvodeMem, maxStep);
00208     // Set the rootfinder function if wanted
00209     if (mCheckForRoots)
00210     {
00211 #if CHASTE_SUNDIALS_VERSION >= 20400
00212         CVodeRootInit(mpCvodeMem, 1, CvodeRootAdaptor);
00213 #else
00214         CVodeRootInit(mpCvodeMem, 1, CvodeRootAdaptor, (void*)(&mData));
00215 #endif
00216     }
00217     // Change max steps if wanted
00218     if (mMaxSteps > 0)
00219     {
00220         CVodeSetMaxNumSteps(mpCvodeMem, mMaxSteps);
00221     }
00222     // Attach a linear solver for Newton iteration
00223     CVDense(mpCvodeMem, rInitialY.size());
00224 }
00225 
00226 void CvodeAdaptor::FreeCvodeMemory()
00227 {
00228 //    assert(!NV_OWN_DATA_STL(mInitialValues));
00229 //    std::vector<double>* pVec = NV_DATA_STL(mInitialValues);
00230 //    double val = (*pVec)[0];
00231     N_VDestroy_Serial(mInitialValues); mInitialValues = NULL;
00232 //    std::cout << "  a: " << val << ", b: " << (*pVec)[0] << std::endl;
00233 
00234     CVodeFree(&mpCvodeMem);
00235 }
00236 
00237 
00238 void CvodeAdaptor::CvodeError(int flag, const char * msg)
00239 {
00240     std::stringstream err;
00241     char* p_flag_name = CVodeGetReturnFlagName(flag);
00242     err << msg << ": " << p_flag_name;
00243     free(p_flag_name);
00244     std::cerr << err.str() << std::endl << std::flush;
00245     EXCEPTION(err.str());
00246 }
00247 
00248 
00249 OdeSolution CvodeAdaptor::Solve(AbstractOdeSystem* pOdeSystem,
00250                                 std::vector<double>& rYValues,
00251                                 double startTime,
00252                                 double endTime,
00253                                 double maxStep,
00254                                 double timeSampling)
00255 {
00256     assert(endTime > startTime);
00257     assert(timeSampling > 0.0);
00258 
00259     mStoppingEventOccurred = false;
00260     if (mCheckForRoots && pOdeSystem->CalculateStoppingEvent(startTime, rYValues) == true)
00261     {
00262         EXCEPTION("(Solve with sampling) Stopping event is true for initial condition");
00263     }
00264 
00265     SetupCvode(pOdeSystem, rYValues, startTime, maxStep);
00266 
00267     TimeStepper stepper(startTime, endTime, timeSampling);
00268     N_Vector yout = mInitialValues;
00269 
00270     // Set up ODE solution
00271     OdeSolution solutions;
00272     solutions.SetNumberOfTimeSteps(stepper.EstimateTimeSteps());
00273     solutions.rGetSolutions().push_back(rYValues);
00274     solutions.rGetTimes().push_back(startTime);
00275     solutions.SetOdeSystemInformation(pOdeSystem->GetSystemInformation());
00276 
00277     // Main time sampling loop
00278     while (!stepper.IsTimeAtEnd() && !mStoppingEventOccurred)
00279     {
00280         double tend;
00281         int ierr = CVode(mpCvodeMem, stepper.GetNextTime(), yout, &tend, CV_NORMAL);
00282         if (ierr<0)
00283         {
00284             FreeCvodeMemory();
00285             CvodeError(ierr, "CVODE failed to solve system");
00286         }
00287         // Store solution
00288         solutions.rGetSolutions().push_back(rYValues);
00289         solutions.rGetTimes().push_back(tend);
00290         if (ierr == CV_ROOT_RETURN)
00291         {
00292             // Stopping event occurred
00293             mStoppingEventOccurred = true;
00294             mStoppingTime = tend;
00295         }
00296         stepper.AdvanceOneTimeStep();
00297     }
00298 
00299     // stepper.EstimateTimeSteps may have been an overestimate...
00300     solutions.SetNumberOfTimeSteps(stepper.GetTotalTimeStepsTaken());
00301 
00302     int ierr = CVodeGetLastStep(mpCvodeMem, &mLastInternalStepSize);
00303     assert(ierr == CV_SUCCESS); ierr=ierr; // avoid unused var warning
00304     FreeCvodeMemory();
00305 
00306     return solutions;
00307 }
00308 
00309 
00310 void CvodeAdaptor::Solve(AbstractOdeSystem* pOdeSystem,
00311                          std::vector<double>& rYValues,
00312                          double startTime,
00313                          double endTime,
00314                          double maxStep)
00315 {
00316     assert(endTime > startTime);
00317 
00318     mStoppingEventOccurred = false;
00319     if (mCheckForRoots && pOdeSystem->CalculateStoppingEvent(startTime, rYValues) == true)
00320     {
00321         EXCEPTION("(Solve) Stopping event is true for initial condition");
00322     }
00323 
00324     SetupCvode(pOdeSystem, rYValues, startTime, maxStep);
00325 
00326     N_Vector yout = mInitialValues;
00327     double tend;
00328     int ierr = CVode(mpCvodeMem, endTime, yout, &tend, CV_NORMAL);
00329     if (ierr<0)
00330     {
00331         FreeCvodeMemory();
00332         CvodeError(ierr, "CVODE failed to solve system");
00333     }
00334     if (ierr == CV_ROOT_RETURN)
00335     {
00336         // Stopping event occurred
00337         mStoppingEventOccurred = true;
00338         mStoppingTime = tend;
00339     }
00340     assert(NV_DATA_S(yout) == &(rYValues[0]));
00341     assert(!NV_OWN_DATA_S(yout));
00342 
00343 //    long int steps;
00344 //    CVodeGetNumSteps(mpCvodeMem, &steps);
00345 //    std::cout << " Solved to " << endTime << " in " << steps << " steps.\n";
00346 
00347     ierr = CVodeGetLastStep(mpCvodeMem, &mLastInternalStepSize);
00348     assert(ierr == CV_SUCCESS); ierr=ierr; // avoid unused var warning
00349     FreeCvodeMemory();
00350 }
00351 
00352 CvodeAdaptor::CvodeAdaptor(double relTol, double absTol)
00353     : AbstractIvpOdeSolver(),
00354       mpCvodeMem(NULL), mInitialValues(NULL),
00355       mRelTol(relTol), mAbsTol(absTol),
00356       mLastInternalStepSize(-0.0),
00357       mMaxSteps(0),
00358       mCheckForRoots(false)
00359 {
00360 }
00361 
00362 void CvodeAdaptor::SetTolerances(double relTol, double absTol)
00363 {
00364     mRelTol = relTol;
00365     mAbsTol = absTol;
00366 }
00367 
00368 double CvodeAdaptor::GetRelativeTolerance()
00369 {
00370     return mRelTol;
00371 }
00372 
00373 double CvodeAdaptor::GetAbsoluteTolerance()
00374 {
00375     return mAbsTol;
00376 }
00377 
00378 double CvodeAdaptor::GetLastStepSize()
00379 {
00380     return mLastInternalStepSize;
00381 }
00382 
00383 void CvodeAdaptor::CheckForStoppingEvents()
00384 {
00385     mCheckForRoots = true;
00386 }
00387 
00388 void CvodeAdaptor::SetMaxSteps(long int numSteps)
00389 {
00390     mMaxSteps = numSteps;
00391 }
00392 
00393 long int CvodeAdaptor::GetMaxSteps()
00394 {
00395     return mMaxSteps;
00396 }
00397 
00398 // Serialization for Boost >= 1.36
00399 #include "SerializationExportWrapperForCpp.hpp"
00400 CHASTE_CLASS_EXPORT(CvodeAdaptor)
00401 
00402 #endif // CHASTE_CVODE

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