OdeSolution.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 "OdeSolution.hpp"
00030 
00031 #include <sstream>
00032 
00033 #include "ColumnDataWriter.hpp"
00034 #include "Exception.hpp"
00035 #include "PetscTools.hpp"
00036 #include "VectorHelperFunctions.hpp"
00037 
00038 OdeSolution::OdeSolution()
00039     : mNumberOfTimeSteps(0u),
00040       mpOdeSystemInformation()
00041 {
00042 }
00043 
00044 
00045 unsigned OdeSolution::GetNumberOfTimeSteps() const
00046 {
00047     return mNumberOfTimeSteps;
00048 }
00049 
00050 
00051 void OdeSolution::SetNumberOfTimeSteps(unsigned numTimeSteps)
00052 {
00053     mNumberOfTimeSteps = numTimeSteps;
00054     mTimes.reserve(numTimeSteps+1);
00055     mSolutions.reserve(numTimeSteps);
00056 }
00057 
00058 
00059 void OdeSolution::SetOdeSystemInformation(boost::shared_ptr<const AbstractOdeSystemInformation> pOdeSystemInfo)
00060 {
00061     mpOdeSystemInformation = pOdeSystemInfo;
00062 }
00063 
00064 std::vector<double> OdeSolution::GetVariableAtIndex(unsigned index) const
00065 {
00066     std::vector<double> answer;
00067     answer.reserve(mTimes.size());
00068     double temp_number;
00069     for (unsigned i=0; i< mTimes.size(); ++i)
00070     {
00071         if (index < mSolutions[0].size())
00072         {
00073             temp_number = mSolutions[i][index];
00074         }
00075         else
00076         {
00077             unsigned offset = mSolutions[0].size();
00078             if (index - offset < mParameters.size())
00079             {
00080                 temp_number = mParameters[index - offset];
00081             }
00082             else
00083             {
00084                 offset += mParameters.size();
00085                 if (index - offset < mDerivedQuantities[0].size())
00086                 {
00087                     temp_number = mDerivedQuantities[i][index - offset];
00088                 }
00089                 else
00090                 {
00091                     EXCEPTION("Invalid index passed to ""GetVariableAtIndex()"".");
00092                 }
00093             }
00094         }
00095         answer.push_back(temp_number);
00096     }
00097     return answer;
00098 }
00099 
00100 std::vector<double> OdeSolution::GetAnyVariable(const std::string& rName) const
00101 {
00102     return GetVariableAtIndex(mpOdeSystemInformation->GetAnyVariableIndex(rName));
00103 }
00104 
00105 std::vector<double>& OdeSolution::rGetTimes()
00106 {
00107     return mTimes;
00108 }
00109 
00110 const std::vector<double>& OdeSolution::rGetTimes() const
00111 {
00112     return mTimes;
00113 }
00114 
00115 std::vector<std::vector<double> >& OdeSolution::rGetSolutions()
00116 {
00117     return mSolutions;
00118 }
00119 
00120 const std::vector<std::vector<double> >& OdeSolution::rGetSolutions() const
00121 {
00122     return mSolutions;
00123 }
00124 
00125 
00126 template<typename VECTOR>
00127 void OdeSolution::CalculateDerivedQuantitiesAndParameters(AbstractParameterisedSystem<VECTOR>* pOdeSystem)
00128 {
00129     assert(pOdeSystem->GetSystemInformation() == mpOdeSystemInformation); // Just in case...
00130     rGetParameters(pOdeSystem);
00131     rGetDerivedQuantities(pOdeSystem);
00132 }
00133 
00134 
00135 template<typename VECTOR>
00136 std::vector<double>& OdeSolution::rGetParameters(AbstractParameterisedSystem<VECTOR>* pOdeSystem)
00137 {
00138     mParameters.clear();
00139     const unsigned num_params = pOdeSystem->GetNumberOfParameters();
00140     if (num_params > 0)
00141     {
00142         mParameters.reserve(num_params);
00143         for (unsigned i=0; i<num_params; ++i)
00144         {
00145             mParameters.push_back(pOdeSystem->GetParameter(i));
00146         }
00147     }
00148     return mParameters;
00149 }
00150 
00151 
00152 std::vector<std::vector<double> >& OdeSolution::rGetDerivedQuantities(AbstractParameterisedSystem<std::vector<double> >* pOdeSystem)
00153 {
00154     assert(pOdeSystem != NULL);
00155     if (mDerivedQuantities.empty() && pOdeSystem->GetNumberOfDerivedQuantities() > 0)
00156     {
00157         assert(mTimes.size() == mSolutions.size()); // Paranoia
00158         mDerivedQuantities.reserve(mTimes.size());
00159         for (unsigned i=0; i<mTimes.size(); i++)
00160         {
00161             mDerivedQuantities.push_back(pOdeSystem->ComputeDerivedQuantities(mTimes[i], mSolutions[i]));
00162         }
00163     }
00164     return mDerivedQuantities;
00165 }
00166 
00167 #ifdef CHASTE_CVODE
00168 std::vector<std::vector<double> >& OdeSolution::rGetDerivedQuantities(AbstractParameterisedSystem<N_Vector>* pOdeSystem)
00169 {
00170     assert(pOdeSystem != NULL);
00171     if (mDerivedQuantities.empty() && pOdeSystem->GetNumberOfDerivedQuantities() > 0)
00172     {
00173         const unsigned num_solutions = mSolutions.size();
00174         assert(mTimes.size() == num_solutions); // Paranoia
00175         mDerivedQuantities.resize(mTimes.size());
00176         N_Vector state_vars = num_solutions > 0 ? N_VNew_Serial(mSolutions[0].size()) : NULL;
00177         for (unsigned i=0; i<num_solutions; i++)
00178         {
00179             CopyFromStdVector(mSolutions[i], state_vars);
00180             N_Vector dqs = pOdeSystem->ComputeDerivedQuantities(mTimes[i], state_vars);
00181             CopyToStdVector(dqs, mDerivedQuantities[i]);
00182             DeleteVector(dqs);
00183         }
00184         DeleteVector(state_vars);
00185     }
00186     assert(mDerivedQuantities.size()==mTimes.size());
00187     return mDerivedQuantities;
00188 }
00189 #endif // CHASTE_CVODE
00190 
00191 
00192 void OdeSolution::WriteToFile(std::string directoryName,
00193                               std::string baseResultsFilename,
00194                               std::string timeUnits,
00195                               unsigned stepsPerRow,
00196                               bool cleanDirectory,
00197                               unsigned precision,
00198                               bool includeDerivedQuantities)
00199 {
00200     assert(stepsPerRow > 0);
00201     assert(mTimes.size() > 0);
00202     assert(mTimes.size() == mSolutions.size());
00203     assert(mpOdeSystemInformation.get() != NULL);
00204     if (mpOdeSystemInformation->GetNumberOfParameters()==0 && mpOdeSystemInformation->GetNumberOfDerivedQuantities() == 0)
00205     {
00206         includeDerivedQuantities = false;
00207     }
00208 
00209     if (includeDerivedQuantities)
00210     {
00211         if ((mDerivedQuantities.empty() || mDerivedQuantities.size()!=mTimes.size()) && mParameters.empty())
00212         {
00213             EXCEPTION("You must first call ""CalculateDerivedQuantitiesAndParameters()"" in order to write derived quantities.");
00214         }
00215     }
00216 
00217     // Write data to a file using ColumnDataWriter
00218     ColumnDataWriter writer(directoryName, baseResultsFilename, cleanDirectory, precision);
00219 
00220     if (!PetscTools::AmMaster())
00221     {
00222         //Only the master actually writes to file
00223         return;
00224     }
00225 
00226     int time_var_id = writer.DefineUnlimitedDimension("Time", timeUnits);
00227 
00228     // Either: the ODE system should have no names&units defined, or it should
00229     // the same number as the number of solutions per timestep.
00230     assert(  mpOdeSystemInformation->rGetStateVariableNames().size()==0 ||
00231             (mpOdeSystemInformation->rGetStateVariableNames().size()==mSolutions[0].size()) );
00232 
00233     unsigned num_vars = mSolutions[0].size();
00234     unsigned num_params = mpOdeSystemInformation->GetNumberOfParameters();
00235     unsigned num_derived_quantities = mpOdeSystemInformation->GetNumberOfDerivedQuantities();
00236 
00237     std::vector<int> var_ids;
00238     var_ids.reserve(num_vars);
00239     if (mpOdeSystemInformation->rGetStateVariableNames().size() > 0)
00240     {
00241         for (unsigned i=0; i<num_vars; i++)
00242         {
00243             var_ids.push_back(writer.DefineVariable(mpOdeSystemInformation->rGetStateVariableNames()[i],
00244                                                     mpOdeSystemInformation->rGetStateVariableUnits()[i]));
00245         }
00246     }
00247     else
00248     {
00249         for (unsigned i=0; i<num_vars; i++)
00250         {
00251             std::stringstream string_stream;
00252             string_stream << "var_" << i;
00253             var_ids.push_back(writer.DefineVariable(string_stream.str(), ""));
00254         }
00255     }
00256 
00257     if (includeDerivedQuantities)
00258     {
00259         var_ids.reserve(num_vars + num_params + num_derived_quantities);
00260         for (unsigned i=0; i<num_params; ++i)
00261         {
00262             var_ids.push_back(writer.DefineVariable(mpOdeSystemInformation->rGetParameterNames()[i],
00263                                                     mpOdeSystemInformation->rGetParameterUnits()[i]));
00264         }
00265         for (unsigned i=0; i<num_derived_quantities; i++)
00266         {
00267             var_ids.push_back(writer.DefineVariable(mpOdeSystemInformation->rGetDerivedQuantityNames()[i],
00268                                                     mpOdeSystemInformation->rGetDerivedQuantityUnits()[i]));
00269         }
00270     }
00271 
00272     if(mSolverName!="")
00273     {
00274         writer.SetCommentForInfoFile("ODE SOLVER: " + mSolverName);
00275     }
00276 
00277     writer.EndDefineMode();
00278 
00279     for (unsigned i=0; i<mSolutions.size(); i+=stepsPerRow)
00280     {
00281         writer.PutVariable(time_var_id, mTimes[i]);
00282         for (unsigned j=0; j<num_vars; j++)
00283         {
00284             writer.PutVariable(var_ids[j], mSolutions[i][j]);
00285         }
00286         if (includeDerivedQuantities)
00287         {
00288             for (unsigned j=0; j<num_params; ++j)
00289             {
00290                 writer.PutVariable(var_ids[j+num_vars], mParameters[j]);
00291             }
00292             for (unsigned j=0; j<num_derived_quantities; j++)
00293             {
00294                 writer.PutVariable(var_ids[j+num_params+num_vars], mDerivedQuantities[i][j]);
00295             }
00296         }
00297         writer.AdvanceAlongUnlimitedDimension();
00298     }
00299     writer.Close();
00300 }
00301 
00302 //
00303 // Explicit instantiation
00304 //
00305 
00310 template std::vector<double>& OdeSolution::rGetParameters(AbstractParameterisedSystem<std::vector<double> >* pOdeSystem);
00311 template void OdeSolution::CalculateDerivedQuantitiesAndParameters(AbstractParameterisedSystem<std::vector<double> >* pOdeSystem);
00312 
00313 #ifdef CHASTE_CVODE
00314 template std::vector<double>& OdeSolution::rGetParameters(AbstractParameterisedSystem<N_Vector>* pOdeSystem);
00315 template void OdeSolution::CalculateDerivedQuantitiesAndParameters(AbstractParameterisedSystem<N_Vector>* pOdeSystem);
00316 #endif // CHASTE_CVODE
00317 

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