Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
AbstractCardiacProblem.hpp
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
37#ifndef ABSTRACTCARDIACPROBLEM_HPP_
38#define ABSTRACTCARDIACPROBLEM_HPP_
39
40#include <string>
41#include <vector>
42#include <cassert>
43#include <climits>
44#include <boost/shared_ptr.hpp>
45
47#include <boost/serialization/vector.hpp>
48#include <boost/serialization/string.hpp>
49#include <boost/serialization/split_member.hpp>
50#include <boost/serialization/shared_ptr.hpp>
51#include "ClassIsAbstract.hpp"
53
54#include "AbstractTetrahedralMesh.hpp"
55#include "AbstractCardiacCell.hpp"
56#include "AbstractCardiacCellFactory.hpp"
57#include "AbstractCardiacTissue.hpp"
58#include "AbstractDynamicLinearPdeSolver.hpp"
59#include "BoundaryConditionsContainer.hpp"
60#include "DistributedVectorFactory.hpp"
61#include "Hdf5DataReader.hpp"
62#include "Hdf5DataWriter.hpp"
63#include "Warnings.hpp"
64#include "AbstractOutputModifier.hpp"
65/*
66 * Archiving extravaganza:
67 *
68 * We archive mesh and tissue through a pointer to an abstract class. All the potential concrete
69 * classes need to be included here, so they are registered with boost. If not, boost won't be
70 * able to find the archiving methods of the concrete class and will throw the following
71 * exception:
72 *
73 * terminate called after throwing an instance of 'boost::archive::archive_exception'
74 * what(): unregistered class
75 *
76 * No member variable is defined to be of any of these clases, so removing them won't
77 * produce any compiler error. The exception above will occur at runtime.
78 *
79 * This might not be even necessary in certain cases, if the file is included implicitly by another
80 * header file or by the test itself. It's safer though.
81 */
82#include "DistributedTetrahedralMesh.hpp"
83#include "TetrahedralMesh.hpp" //May be needed for unarchiving a mesh
84#include "MonodomainTissue.hpp"
85#include "BidomainTissue.hpp"
86
90class AbstractUntemplatedCardiacProblem : private boost::noncopyable
91{
92public:
96};
97
112template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
114{
115 friend class TestBidomainWithBath;
116 friend class TestMonodomainProblem;
117 friend class TestCardiacSimulationArchiver;
118
120 typedef typename boost::shared_ptr<BoundaryConditionsContainer<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM> >
122
123private:
126
133 template<class Archive>
134 void save(Archive & archive, const unsigned int version) const
135 {
136 if (version >= 1)
137 {
138 const unsigned element_dim=ELEMENT_DIM;
139 archive & element_dim;
140 const unsigned space_dim=SPACE_DIM;
141 archive & space_dim;
142 const unsigned problem_dim=PROBLEM_DIM;
143 archive & problem_dim;
144 }
145 archive & mMeshFilename;
146 archive & mpMesh;
147 //archive & mAllocatedMemoryForMesh; // Mesh is deleted by AbstractCardiacTissue
148
149 // We shouldn't ever have to save the old version
150 assert(version >= 2);
151// {
152// bool use_matrix_based_assembly = true;
153// archive & use_matrix_based_assembly;
154// }
155
156 archive & mWriteInfo;
157 archive & mPrintOutput;
158 archive & mNodesToOutput;
159 //archive & mVoltageColumnId; // Created by InitialiseWriter, called from Solve
160 //archive & mExtraVariablesId; // Created by InitialiseWriter, called from Solve
161 //archive & mTimeColumnId; // Created by InitialiseWriter, called from Solve
162 //archive & mNodeColumnId; // Created by InitialiseWriter, called from Solve
163 //archive & mpWriter; // Created by InitialiseWriter, called from Solve
164 archive & mpCardiacTissue;
165 //archive & mpSolver; // Only exists during calls to the Solve method
166 bool has_solution = (mSolution != NULL);
167 archive & has_solution;
168 if (has_solution)
169 {
172 Hdf5DataWriter writer(*mpMesh->GetDistributedVectorFactory(), ArchiveLocationInfo::GetArchiveRelativePath(), "AbstractCardiacProblem_mSolution", false);
173 writer.DefineFixedDimension(mpMesh->GetDistributedVectorFactory()->GetProblemSize());
174 writer.DefineUnlimitedDimension("Time", "msec", 1);
175
176 int vm_col = writer.DefineVariable("Vm","mV");
177
178 if (PROBLEM_DIM==1)
179 {
180 writer.EndDefineMode();
181 writer.PutUnlimitedVariable(0.0);
182 writer.PutVector(vm_col, mSolution);
183 }
184
185 if (PROBLEM_DIM==2)
186 {
187 int phie_col = writer.DefineVariable("Phie","mV");
188 std::vector<int> variable_ids;
189 variable_ids.push_back(vm_col);
190 variable_ids.push_back(phie_col);
191 writer.EndDefineMode();
192 writer.PutUnlimitedVariable(0.0);
193 writer.PutStripedVector(variable_ids, mSolution);
194 }
195
196 writer.Close();
197
198 }
199 archive & mCurrentTime;
200
201 // Save boundary conditions
204
205 if (version >= 3)
206 {
207 archive & mOutputModifiers;
208 }
209
210 if (version >= 4)
211 {
212 archive & mUseHdf5DataWriterCache;
214 }
215 }
216
223 template<class Archive>
224 void load(Archive & archive, const unsigned int version)
225 {
226 if (version >= 1)
227 {
228 unsigned element_dim;
229 unsigned space_dim;
230 unsigned problem_dim;
231 archive & element_dim;
232 archive & space_dim;
233 archive & problem_dim;
234 if ((element_dim != ELEMENT_DIM) ||(space_dim != SPACE_DIM) ||(problem_dim != PROBLEM_DIM))
235 {
236 /*If we carry on from this point then the mesh produced by unarchiving from the
237 * archive is templated as AbstractTetrahedralMesh<element_dim, space_dim>
238 * which doesn't match AbstractTetrahedralMesh<ELEMENT_DIM, SPACE_DIM>* mpMesh.
239 * Boost will through away the unarchived one, without deleting it properly and
240 * then set mpMesh=NULL. We need to avoid this happening by bailing out.
241 */
242 EXCEPTION("Failed to load from checkpoint because the dimensions of the archive do not match the object it's being read into.");
243 }
244 }
245 archive & mMeshFilename;
246 archive & mpMesh;
247 assert(mpMesh != NULL); //If NULL then loading mesh has failed without an exception so Boost has given up on the mesh. This would happen if a 2-dimensional mesh was successfully unarchived but mpMesh was expecting a 3-d mesh etc.
248 //archive & mAllocatedMemoryForMesh; // Will always be true after a load
249
250 if (version < 2)
251 {
252 bool use_matrix_based_assembly;
253 archive & use_matrix_based_assembly;
254 }
255
256 archive & mWriteInfo;
257 archive & mPrintOutput;
258 archive & mNodesToOutput;
259 //archive & mVoltageColumnId; // Created by InitialiseWriter, called from Solve
260 //archive & mExtraVariablesId; // Created by InitialiseWriter, called from Solve
261 //archive & mTimeColumnId; // Created by InitialiseWriter, called from Solve
262 //archive & mNodeColumnId; // Created by InitialiseWriter, called from Solve
263 //archive & mpWriter; // Created by InitialiseWriter, called from Solve
264 archive & mpCardiacTissue;
265 //archive & mpSolver; // Only exists during calls to the Solve method
266 bool has_solution;
267 archive & has_solution;
268 if ((has_solution) && PROBLEM_DIM < 3)
269 {
272 mSolution = mpMesh->GetDistributedVectorFactory()->CreateVec(PROBLEM_DIM);
273 DistributedVector mSolution_distri = mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(mSolution);
274
275 Vec vm = mpMesh->GetDistributedVectorFactory()->CreateVec();
276 Vec phie = mpMesh->GetDistributedVectorFactory()->CreateVec();
277
278 std::string archive_dir = ArchiveLocationInfo::GetArchiveRelativePath();
279 Hdf5DataReader reader(archive_dir, "AbstractCardiacProblem_mSolution", !FileFinder::IsAbsolutePath(archive_dir));
280 reader.GetVariableOverNodes(vm, "Vm", 0);
281
282 if (PROBLEM_DIM==1)
283 {
284 //reader.Close(); // no need to call close explicitly, done in the destructor
285
286 DistributedVector vm_distri = mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(vm);
287
288 DistributedVector::Stripe mSolution_vm(mSolution_distri,0);
289
290 for (DistributedVector::Iterator index = mSolution_distri.Begin();
291 index != mSolution_distri.End();
292 ++index)
293 {
294 mSolution_vm[index] = vm_distri[index];
295 }
296 }
297
298 if (PROBLEM_DIM==2)
299 {
300 reader.GetVariableOverNodes(phie, "Phie", 0);
301 //reader.Close(); // no need to call close explicitly, done in the destructor
302
303 DistributedVector vm_distri = mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(vm);
304 DistributedVector phie_distri = mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(phie);
305
306 DistributedVector::Stripe mSolution_vm(mSolution_distri,0);
307 DistributedVector::Stripe mSolution_phie(mSolution_distri,1);
308
309 for (DistributedVector::Iterator index = mSolution_distri.Begin();
310 index != mSolution_distri.End();
311 ++index)
312 {
313 mSolution_vm[index] = vm_distri[index];
314 mSolution_phie[index] = phie_distri[index];
315 }
316 }
317
318 mSolution_distri.Restore();
319
322
323 }
324 archive & mCurrentTime;
325
326 // Load boundary conditions
329
330 if (version >= 3)
331 {
332 archive & mOutputModifiers;
333 }
334
335 if (version >= 4)
336 {
337 archive & mUseHdf5DataWriterCache;
339 }
340 }
341
342 BOOST_SERIALIZATION_SPLIT_MEMBER()
343
344
351 template<class Archive>
352 void SaveBoundaryConditions(Archive & archive,
353 AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
354 BccType pBcc) const
355 {
357 }
358
366 template<class Archive>
368 Archive & archive,
370 {
371 // Load pointer from archive
372 BccType p_bcc;
374
375 // Fill in the conditions, if we have a container and it's not already full
376 if (p_bcc)
377 {
378 p_bcc->LoadFromArchive(*ProcessSpecificArchive<Archive>::Get(), pMesh);
379 }
380
381 return p_bcc;
382 }
383
384protected:
387 std::string mMeshFilename;
388
395
397 std::vector<unsigned> mNodesToOutput;
398
402 std::vector<unsigned> mExtraVariablesId;
407
410
421
425
434
437
445
453
461 virtual void CreateMeshFromHeartConfig();
462
466 template<unsigned DIM, unsigned ELEC_PROB_DIM>
468
473
478
483
487 std::vector<boost::shared_ptr<AbstractOutputModifier> > mOutputModifiers;
488
489public:
496
501
505 virtual ~AbstractCardiacProblem();
506
514 void Initialise();
515
521 void SetNodesPerProcessorFilename(const std::string& rFilename);
522
528
536 virtual void PreSolveChecks();
537
551 virtual Vec CreateInitialCondition();
552
559
565 void PrintOutput(bool rPrintOutput);
566
572 void SetWriteInfo(bool writeInfo = true);
573
585
592
596 double GetCurrentTime();
597
602
607
617 void Solve();
618
627
635 virtual void WriteInfo(double time)=0;
636
641 virtual void DefineWriterColumns(bool extending);
642
647 void DefineExtraVariablesWriterColumns(bool extending);
648
655 virtual void WriteOneStep(double time, Vec voltageVec) = 0;
656
661
672 bool InitialiseWriter();
673
680 void SetUseHdf5DataWriterCache(bool useCache=true);
681
702
712 void SetOutputNodes(std::vector<unsigned> & rNodesToOutput);
713
718
726 virtual void AtBeginningOfTimestep(double time)
727 {}
728
736 virtual void OnEndOfTimestep(double time)
737 {}
738
746 virtual void SetUpAdditionalStoppingTimes(std::vector<double>& rAdditionalStoppingTimes)
747 {}
748
750 // and no controller, in which case the default is used.
751
757 void SetUseTimeAdaptivityController(bool useAdaptivity,
758 AbstractTimeAdaptivityController* pController = NULL);
759
780 template<class Archive>
781 void LoadExtraArchive(Archive & archive, unsigned version);
782
786 virtual bool GetHasBath();
787
792 virtual void SetElectrodes();
793
799 void AddOutputModifier( boost::shared_ptr<AbstractOutputModifier> pOutputModifier)
800 {
801 mOutputModifiers.push_back(pOutputModifier);
802 }
803};
804
806
807
808template<unsigned DIM>
809class BidomainProblem;
810
811template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
812template<class Archive>
814{
815 // The vector factory must be loaded, but isn't needed for anything.
816 DistributedVectorFactory* p_mesh_factory;
817 archive >> p_mesh_factory;
818
819 // How many processes were used by the saving simulation?
820 DistributedVectorFactory* p_original_factory = p_mesh_factory->GetOriginalFactory();
821 unsigned orig_num_procs = 1;
822 if (p_original_factory)
823 {
824 orig_num_procs = p_original_factory->GetNumProcs();
825 }
826
827 // The cardiac cells - load only the cells we actually own
828 mpCardiacTissue->LoadCardiacCells(archive, version);
829
830 {
831 DistributedVectorFactory* p_pde_factory;
832 archive >> p_pde_factory;
833 assert(p_pde_factory == p_mesh_factory); // Paranoia...
834 }
835 // We no longer need this vector factory, since we already have our own.
836 delete p_mesh_factory;
837
838 // The boundary conditions
839 BccType p_bcc;
840 archive >> p_bcc;
841 if (p_bcc)
842 {
844 {
846 mpBoundaryConditionsContainer->LoadFromArchive(archive, mpMesh);
847 }
848 else
849 {
850 // If the mesh which was archived was a TetrahedralMesh then we have all the boundary conditions
851 // in every process-specific archive. We no longer test for this.
852// LCOV_EXCL_START
853 if (!dynamic_cast<DistributedTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>*>(mpMesh) && orig_num_procs > 1)
854 {
855 // The correct way to do this should be:
856 // p_bcc->LoadFromArchive(archive, mpMesh);
857 WARNING("Loading from a parallel archive which used a non-distributed mesh. This scenario should work but is not fully tested.");
858 }
859// LCOV_EXCL_STOP
860 mpBoundaryConditionsContainer->MergeFromArchive(archive, mpMesh);
861 }
862 }
863 BccType p_default_bcc;
864 archive >> p_default_bcc;
865 if (p_default_bcc)
866 {
867 // This always holds, so we never need to load the default BCs
868 assert(p_bcc == p_default_bcc);
869 }
870
871 // Are we a bidomain problem?
872 BidomainProblem<ELEMENT_DIM>* p_bidomain_problem = dynamic_cast<BidomainProblem<ELEMENT_DIM>*>(this);
873 if (p_bidomain_problem)
874 {
875 assert(ELEMENT_DIM == SPACE_DIM);
876 p_bidomain_problem->LoadExtraArchiveForBidomain(archive, version);
877 }
878}
879
880namespace boost {
881namespace serialization {
888template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
889struct version<AbstractCardiacProblem<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM> >
890{
893};
894} // namespace serialization
895} // namespace boost
896#endif /*ABSTRACTCARDIACPROBLEM_HPP_*/
gcov doesn't like this file...
#define TEMPLATED_CLASS_IS_ABSTRACT_3_UNSIGNED(T)
#define EXCEPTION(message)
AbstractCardiacTissue< ELEMENT_DIM, SPACE_DIM > * mpCardiacTissue
void SetWriteInfo(bool writeInfo=true)
virtual AbstractDynamicLinearPdeSolver< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > * CreateSolver()=0
virtual void OnEndOfTimestep(double time)
void load(Archive &archive, const unsigned int version)
std::vector< boost::shared_ptr< AbstractOutputModifier > > mOutputModifiers
virtual AbstractCardiacTissue< ELEMENT_DIM, SPACE_DIM > * CreateCardiacTissue()=0
void AddOutputModifier(boost::shared_ptr< AbstractOutputModifier > pOutputModifier)
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
BccType LoadBoundaryConditions(Archive &archive, AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh)
void SetOutputNodes(std::vector< unsigned > &rNodesToOutput)
void DefineExtraVariablesWriterColumns(bool extending)
DistributedVector GetSolutionDistributedVector()
void save(Archive &archive, const unsigned int version) const
void SetHdf5DataWriterTargetChunkSizeAndAlignment(hsize_t size)
virtual void SetUpAdditionalStoppingTimes(std::vector< double > &rAdditionalStoppingTimes)
AbstractTimeAdaptivityController * mpTimeAdaptivityController
void SetUseTimeAdaptivityController(bool useAdaptivity, AbstractTimeAdaptivityController *pController=NULL)
void SetMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh)
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > & rGetMesh()
AbstractDynamicLinearPdeSolver< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > * mpSolver
boost::shared_ptr< BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > > BccType
void SaveBoundaryConditions(Archive &archive, AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, BccType pBcc) const
virtual void WriteOneStep(double time, Vec voltageVec)=0
virtual void WriteInfo(double time)=0
void PrintOutput(bool rPrintOutput)
std::vector< unsigned > mExtraVariablesId
void SetBoundaryConditionsContainer(BccType pBcc)
friend class boost::serialization::access
AbstractCardiacCellFactory< ELEMENT_DIM, SPACE_DIM > * mpCellFactory
std::vector< unsigned > mNodesToOutput
AbstractCardiacTissue< ELEMENT_DIM, SPACE_DIM > * GetTissue()
void SetUseHdf5DataWriterCache(bool useCache=true)
void LoadExtraArchive(Archive &archive, unsigned version)
void SetNodesPerProcessorFilename(const std::string &rFilename)
virtual void DefineWriterColumns(bool extending)
virtual void AtBeginningOfTimestep(double time)
static std::string GetArchiveRelativePath()
void LoadExtraArchiveForBidomain(Archive &archive, unsigned version)
DistributedVectorFactory * GetOriginalFactory()
static bool IsAbsolutePath(const std::string &rPath)
void GetVariableOverNodes(Vec data, const std::string &rVariableName, unsigned timestep=0)
void DefineUnlimitedDimension(const std::string &rVariableName, const std::string &rVariableUnits, unsigned estimatedLength=1)
void PutUnlimitedVariable(double value)
void DefineFixedDimension(long dimensionSize)
int DefineVariable(const std::string &rVariableName, const std::string &rVariableUnits)
virtual void EndDefineMode()
void PutVector(int variableID, Vec petscVector)
void PutStripedVector(std::vector< int > variableIDs, Vec petscVector)
static void Destroy(Vec &rVec)
CHASTE_VERSION_CONTENT(4)
Macro to set the version number of templated archive in known versions of Boost.