Chaste  Release::2018.1
OdeSolution.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "OdeSolution.hpp"
37 
38 #include <sstream>
39 
40 #include "ColumnDataWriter.hpp"
41 #include "Exception.hpp"
42 #include "PetscTools.hpp"
44 
46  : mNumberOfTimeSteps(0u),
47  mpOdeSystemInformation()
48 {
49 }
50 
51 
53 {
54  return mNumberOfTimeSteps;
55 }
56 
57 
58 void OdeSolution::SetNumberOfTimeSteps(unsigned numTimeSteps)
59 {
60  mNumberOfTimeSteps = numTimeSteps;
61  mTimes.reserve(numTimeSteps+1);
62  mSolutions.reserve(numTimeSteps);
63 }
64 
65 
66 void OdeSolution::SetOdeSystemInformation(boost::shared_ptr<const AbstractOdeSystemInformation> pOdeSystemInfo)
67 {
68  mpOdeSystemInformation = pOdeSystemInfo;
69 }
70 
71 std::vector<double> OdeSolution::GetVariableAtIndex(unsigned index) const
72 {
73  std::vector<double> answer;
74  answer.reserve(mTimes.size());
75  double temp_number;
76  for (unsigned i=0; i< mTimes.size(); ++i)
77  {
78  if (index < mSolutions[0].size())
79  {
80  temp_number = mSolutions[i][index];
81  }
82  else
83  {
84  unsigned offset = mSolutions[0].size();
85  if (index - offset < mParameters.size())
86  {
87  temp_number = mParameters[index - offset];
88  }
89  else
90  {
91  offset += mParameters.size();
92  if (index - offset < mDerivedQuantities[0].size())
93  {
94  temp_number = mDerivedQuantities[i][index - offset];
95  }
96  else
97  {
98  EXCEPTION("Invalid index passed to ""GetVariableAtIndex()"".");
99  }
100  }
101  }
102  answer.push_back(temp_number);
103  }
104  return answer;
105 }
106 
107 std::vector<double> OdeSolution::GetAnyVariable(const std::string& rName) const
108 {
109  return GetVariableAtIndex(mpOdeSystemInformation->GetAnyVariableIndex(rName));
110 }
111 
112 std::vector<double>& OdeSolution::rGetTimes()
113 {
114  return mTimes;
115 }
116 
117 const std::vector<double>& OdeSolution::rGetTimes() const
118 {
119  return mTimes;
120 }
121 
122 std::vector<std::vector<double> >& OdeSolution::rGetSolutions()
123 {
124  return mSolutions;
125 }
126 
127 const std::vector<std::vector<double> >& OdeSolution::rGetSolutions() const
128 {
129  return mSolutions;
130 }
131 
132 
133 template<typename VECTOR>
135 {
136  assert(pOdeSystem->GetSystemInformation() == mpOdeSystemInformation); // Just in case...
137  rGetParameters(pOdeSystem);
138  rGetDerivedQuantities(pOdeSystem);
139 }
140 
141 
142 template<typename VECTOR>
144 {
145  mParameters.clear();
146  const unsigned num_params = pOdeSystem->GetNumberOfParameters();
147  if (num_params > 0)
148  {
149  mParameters.reserve(num_params);
150  for (unsigned i=0; i<num_params; ++i)
151  {
152  mParameters.push_back(pOdeSystem->GetParameter(i));
153  }
154  }
155  return mParameters;
156 }
157 
158 
159 std::vector<std::vector<double> >& OdeSolution::rGetDerivedQuantities(AbstractParameterisedSystem<std::vector<double> >* pOdeSystem)
160 {
161  assert(pOdeSystem != nullptr);
162  if (mDerivedQuantities.empty() && pOdeSystem->GetNumberOfDerivedQuantities() > 0)
163  {
164  assert(mTimes.size() == mSolutions.size()); // Paranoia
165  mDerivedQuantities.reserve(mTimes.size());
166  for (unsigned i=0; i<mTimes.size(); i++)
167  {
168  mDerivedQuantities.push_back(pOdeSystem->ComputeDerivedQuantities(mTimes[i], mSolutions[i]));
169  }
170  }
171  return mDerivedQuantities;
172 }
173 
174 #ifdef CHASTE_CVODE
175 std::vector<std::vector<double> >& OdeSolution::rGetDerivedQuantities(AbstractParameterisedSystem<N_Vector>* pOdeSystem)
176 {
177  assert(pOdeSystem != nullptr);
178  if (mDerivedQuantities.empty() && pOdeSystem->GetNumberOfDerivedQuantities() > 0)
179  {
180  const unsigned num_solutions = mSolutions.size();
181  assert(mTimes.size() == num_solutions); // Paranoia
182  mDerivedQuantities.resize(mTimes.size());
183  N_Vector state_vars = num_solutions > 0 ? N_VNew_Serial(mSolutions[0].size()) : nullptr;
184  for (unsigned i=0; i<num_solutions; i++)
185  {
186  CopyFromStdVector(mSolutions[i], state_vars);
187  N_Vector dqs = pOdeSystem->ComputeDerivedQuantities(mTimes[i], state_vars);
189  DeleteVector(dqs);
190  }
191  DeleteVector(state_vars);
192  }
193  assert(mDerivedQuantities.size()==mTimes.size());
194  return mDerivedQuantities;
195 }
196 #endif // CHASTE_CVODE
197 
198 
199 void OdeSolution::WriteToFile(std::string directoryName,
200  std::string baseResultsFilename,
201  std::string timeUnits,
202  unsigned stepsPerRow,
203  bool cleanDirectory,
204  unsigned precision,
205  bool includeDerivedQuantities)
206 {
207  assert(stepsPerRow > 0);
208  assert(mTimes.size() > 0);
209  assert(mTimes.size() == mSolutions.size());
210  assert(mpOdeSystemInformation.get() != nullptr);
211  if (mpOdeSystemInformation->GetNumberOfParameters()==0 && mpOdeSystemInformation->GetNumberOfDerivedQuantities() == 0)
212  {
213  includeDerivedQuantities = false;
214  }
215 
216  if (includeDerivedQuantities)
217  {
218  if ((mDerivedQuantities.empty() || mDerivedQuantities.size()!=mTimes.size()) && mParameters.empty())
219  {
220  EXCEPTION("You must first call ""CalculateDerivedQuantitiesAndParameters()"" in order to write derived quantities.");
221  }
222  }
223 
224  // Write data to a file using ColumnDataWriter
225  ColumnDataWriter writer(directoryName, baseResultsFilename, cleanDirectory, precision);
226 
227  if (!PetscTools::AmMaster())
228  {
229  //Only the master actually writes to file
230  return;
231  }
232 
233  int time_var_id = writer.DefineUnlimitedDimension("Time", timeUnits);
234 
235  // Either: the ODE system should have no names&units defined, or it should
236  // the same number as the number of solutions per timestep.
237  assert( mpOdeSystemInformation->rGetStateVariableNames().size()==0 ||
238  (mpOdeSystemInformation->rGetStateVariableNames().size()==mSolutions[0].size()) );
239 
240  unsigned num_vars = mSolutions[0].size();
241  unsigned num_params = mpOdeSystemInformation->GetNumberOfParameters();
242  unsigned num_derived_quantities = mpOdeSystemInformation->GetNumberOfDerivedQuantities();
243 
244  std::vector<int> var_ids;
245  var_ids.reserve(num_vars);
246  if (mpOdeSystemInformation->rGetStateVariableNames().size() > 0)
247  {
248  for (unsigned i=0; i<num_vars; i++)
249  {
250  var_ids.push_back(writer.DefineVariable(mpOdeSystemInformation->rGetStateVariableNames()[i],
251  mpOdeSystemInformation->rGetStateVariableUnits()[i]));
252  }
253  }
254  else
255  {
256  for (unsigned i=0; i<num_vars; i++)
257  {
258  std::stringstream string_stream;
259  string_stream << "var_" << i;
260  var_ids.push_back(writer.DefineVariable(string_stream.str(), ""));
261  }
262  }
263 
264  if (includeDerivedQuantities)
265  {
266  var_ids.reserve(num_vars + num_params + num_derived_quantities);
267  for (unsigned i=0; i<num_params; ++i)
268  {
269  var_ids.push_back(writer.DefineVariable(mpOdeSystemInformation->rGetParameterNames()[i],
270  mpOdeSystemInformation->rGetParameterUnits()[i]));
271  }
272  for (unsigned i=0; i<num_derived_quantities; i++)
273  {
274  var_ids.push_back(writer.DefineVariable(mpOdeSystemInformation->rGetDerivedQuantityNames()[i],
275  mpOdeSystemInformation->rGetDerivedQuantityUnits()[i]));
276  }
277  }
278 
279  if (mSolverName != "")
280  {
281  writer.SetCommentForInfoFile("ODE SOLVER: " + mSolverName);
282  }
283 
284  writer.EndDefineMode();
285 
286  for (unsigned i=0; i<mSolutions.size(); i+=stepsPerRow)
287  {
288  writer.PutVariable(time_var_id, mTimes[i]);
289  for (unsigned j=0; j<num_vars; j++)
290  {
291  writer.PutVariable(var_ids[j], mSolutions[i][j]);
292  }
293  if (includeDerivedQuantities)
294  {
295  for (unsigned j=0; j<num_params; ++j)
296  {
297  writer.PutVariable(var_ids[j+num_vars], mParameters[j]);
298  }
299  for (unsigned j=0; j<num_derived_quantities; j++)
300  {
301  writer.PutVariable(var_ids[j+num_params+num_vars], mDerivedQuantities[i][j]);
302  }
303  }
305  }
306  writer.Close();
307 }
308 
309 // Explicit instantiation
310 
315 template std::vector<double>& OdeSolution::rGetParameters(AbstractParameterisedSystem<std::vector<double> >* pOdeSystem);
316 template void OdeSolution::CalculateDerivedQuantitiesAndParameters(AbstractParameterisedSystem<std::vector<double> >* pOdeSystem);
317 
318 #ifdef CHASTE_CVODE
319 template std::vector<double>& OdeSolution::rGetParameters(AbstractParameterisedSystem<N_Vector>* pOdeSystem);
321 #endif // CHASTE_CVODE
322 
boost::shared_ptr< const AbstractOdeSystemInformation > mpOdeSystemInformation
Definition: OdeSolution.hpp:84
unsigned mNumberOfTimeSteps
Definition: OdeSolution.hpp:62
unsigned GetNumberOfTimeSteps() const
Definition: OdeSolution.cpp:52
std::vector< std::vector< double > > & rGetSolutions()
virtual void PutVariable(int variableID, double variableValue, long dimensionPosition=-1)
std::vector< std::vector< double > > & rGetDerivedQuantities(AbstractParameterisedSystem< std::vector< double > > *pOdeSystem)
std::vector< std::vector< double > > mDerivedQuantities
Definition: OdeSolution.hpp:71
int DefineUnlimitedDimension(const std::string &rDimensionName, const std::string &rDimensionUnits)
#define EXCEPTION(message)
Definition: Exception.hpp:143
virtual void AdvanceAlongUnlimitedDimension()
static bool AmMaster()
Definition: PetscTools.cpp:120
std::vector< double > mParameters
Definition: OdeSolution.hpp:74
virtual void Close()
virtual void EndDefineMode()
double GetParameter(unsigned index) const
std::vector< double > & rGetParameters(AbstractParameterisedSystem< VECTOR > *pOdeSystem)
void SetOdeSystemInformation(boost::shared_ptr< const AbstractOdeSystemInformation > pOdeSystemInfo)
Definition: OdeSolution.cpp:66
std::vector< double > GetVariableAtIndex(unsigned index) const
Definition: OdeSolution.cpp:71
std::string mSolverName
Definition: OdeSolution.hpp:77
void WriteToFile(std::string directoryName, std::string baseResultsFilename, std::string timeUnits, unsigned stepsPerRow=1, bool cleanDirectory=true, unsigned precision=8, bool includeDerivedQuantities=false)
void CalculateDerivedQuantitiesAndParameters(AbstractParameterisedSystem< VECTOR > *pOdeSystem)
boost::shared_ptr< const AbstractOdeSystemInformation > GetSystemInformation() const
int DefineVariable(const std::string &rVariableName, const std::string &rVariableUnits)
std::vector< std::vector< double > > mSolutions
Definition: OdeSolution.hpp:68
std::vector< double > & rGetTimes()
std::vector< double > mTimes
Definition: OdeSolution.hpp:65
void SetCommentForInfoFile(std::string comment)
virtual VECTOR ComputeDerivedQuantities(double time, const VECTOR &rState)
void CopyFromStdVector(const std::vector< double > &rSrc, VECTOR &rDest)
void CopyToStdVector(const VECTOR &rSrc, std::vector< double > &rDest)
std::vector< double > GetAnyVariable(const std::string &rName) const
void SetNumberOfTimeSteps(unsigned numTimeSteps)
Definition: OdeSolution.cpp:58
void DeleteVector(VECTOR &rVec)