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