Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
AbstractCardiacProblem.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#include "AbstractCardiacProblem.hpp"
36
37#include "DistributedVector.hpp"
38#include "Exception.hpp"
39#include "GenericMeshReader.hpp"
40#include "Hdf5ToCmguiConverter.hpp"
41#include "Hdf5ToMeshalyzerConverter.hpp"
42#include "Hdf5ToVtkConverter.hpp"
43#include "HeartConfig.hpp"
44#include "HeartEventHandler.hpp"
45#include "LinearSystem.hpp"
46#include "PetscTools.hpp"
47#include "PostProcessingWriter.hpp"
48#include "ProgressReporter.hpp"
49#include "TimeStepper.hpp"
50
51template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
54 : mMeshFilename(""), // i.e. undefined
55 mAllocatedMemoryForMesh(false),
56 mWriteInfo(false),
57 mPrintOutput(true),
58 mpCardiacTissue(NULL),
59 mpSolver(NULL),
60 mpCellFactory(pCellFactory),
61 mpMesh(NULL),
62 mSolution(NULL),
63 mCurrentTime(0.0),
64 mpTimeAdaptivityController(NULL),
65 mpWriter(NULL),
66 mUseHdf5DataWriterCache(false),
67 mHdf5DataWriterChunkSizeAndAlignment(0)
68{
69 assert(mNodesToOutput.empty());
70 if (!mpCellFactory)
71 {
72 EXCEPTION("AbstractCardiacProblem: Please supply a cell factory pointer to your cardiac problem constructor.");
73 }
74 HeartEventHandler::BeginEvent(HeartEventHandler::EVERYTHING);
75}
76
77template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
79 // It doesn't really matter what we initialise these to, as they'll be overwritten by
80 // the serialization methods
81 : mMeshFilename(""),
82 mAllocatedMemoryForMesh(false), // Handled by AbstractCardiacTissue
83 mWriteInfo(false),
84 mPrintOutput(true),
85 mVoltageColumnId(UINT_MAX),
86 mTimeColumnId(UINT_MAX),
87 mNodeColumnId(UINT_MAX),
88 mpCardiacTissue(NULL),
89 mpSolver(NULL),
90 mpCellFactory(NULL),
91 mpMesh(NULL),
92 mSolution(NULL),
93 mCurrentTime(0.0),
94 mpTimeAdaptivityController(NULL),
95 mpWriter(NULL),
96 mUseHdf5DataWriterCache(false),
97 mHdf5DataWriterChunkSizeAndAlignment(0)
98{
99}
100
101template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
103{
104 delete mpCardiacTissue;
105 if (mSolution)
106 {
107 PetscTools::Destroy(mSolution);
108 }
109
110 if (mAllocatedMemoryForMesh)
111 {
112 delete mpMesh;
113 }
114}
115
116template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
118{
119 HeartEventHandler::BeginEvent(HeartEventHandler::READ_MESH);
120 if (mpMesh)
121 {
123 {
124 WARNING("Using a non-distributed mesh in a parallel simulation is not a good idea.");
125 }
126 }
127 else
128 {
129 // If no mesh has been passed, we get it from the configuration file
130 try
131 {
132 if (HeartConfig::Instance()->GetLoadMesh())
133 {
134 CreateMeshFromHeartConfig();
135 std::shared_ptr<AbstractMeshReader<ELEMENT_DIM, SPACE_DIM> > p_mesh_reader
136 = GenericMeshReader<ELEMENT_DIM, SPACE_DIM>(HeartConfig::Instance()->GetMeshName());
137 mpMesh->ConstructFromMeshReader(*p_mesh_reader);
138 }
139 else if (HeartConfig::Instance()->GetCreateMesh())
140 {
141 CreateMeshFromHeartConfig();
142 assert(HeartConfig::Instance()->GetSpaceDimension() == SPACE_DIM);
143 double inter_node_space = HeartConfig::Instance()->GetInterNodeSpace();
144
145 switch (HeartConfig::Instance()->GetSpaceDimension())
146 {
147 case 1:
148 {
149 c_vector<double, 1> fibre_length;
150 HeartConfig::Instance()->GetFibreLength(fibre_length);
151 mpMesh->ConstructRegularSlabMesh(inter_node_space, fibre_length[0]);
152 break;
153 }
154 case 2:
155 {
156 c_vector<double, 2> sheet_dimensions; //cm
157 HeartConfig::Instance()->GetSheetDimensions(sheet_dimensions);
158 mpMesh->ConstructRegularSlabMesh(inter_node_space, sheet_dimensions[0], sheet_dimensions[1]);
159 break;
160 }
161 case 3:
162 {
163 c_vector<double, 3> slab_dimensions; //cm
164 HeartConfig::Instance()->GetSlabDimensions(slab_dimensions);
165 mpMesh->ConstructRegularSlabMesh(inter_node_space, slab_dimensions[0], slab_dimensions[1], slab_dimensions[2]);
166 break;
167 }
168 default:
170 }
171 }
172 else
173 {
175 }
176
177 mAllocatedMemoryForMesh = true;
178 }
179 catch (Exception& e)
180 {
181 EXCEPTION(std::string("No mesh given: define it in XML parameters file or call SetMesh()\n") + e.GetShortMessage());
182 }
183 }
184 mpCellFactory->SetMesh(mpMesh);
185 HeartEventHandler::EndEvent(HeartEventHandler::READ_MESH);
186
187 HeartEventHandler::BeginEvent(HeartEventHandler::INITIALISE);
188
189 // If the user requested transmural stuff, we fill in the mCellHeterogeneityAreas here
190 if (HeartConfig::Instance()->AreCellularTransmuralHeterogeneitiesRequested())
191 {
192 mpCellFactory->FillInCellularTransmuralAreas();
193 }
194
195 delete mpCardiacTissue; // In case we're called twice
196 mpCardiacTissue = CreateCardiacTissue();
197
198 HeartEventHandler::EndEvent(HeartEventHandler::INITIALISE);
199
200 // Delete any previous solution, so we get a fresh initial condition
201 if (mSolution)
202 {
203 HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
204 PetscTools::Destroy(mSolution);
205 mSolution = NULL;
206 HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
207 }
208
209 // Always start at time zero
210 mCurrentTime = 0.0;
211
212 // For Bidomain with bath, this is where we set up the electrodes
213 SetElectrodes();
214}
215
216template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
221
222template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
227
228template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
230{
231 if (mpCardiacTissue == NULL) // if tissue is NULL, Initialise() probably hasn't been called
232 {
233 EXCEPTION("Cardiac tissue is null, Initialise() probably hasn't been called");
234 }
235 if (HeartConfig::Instance()->GetSimulationDuration() <= mCurrentTime)
236 {
237 EXCEPTION("End time should be in the future");
238 }
239 if (mPrintOutput)
240 {
241 if ((HeartConfig::Instance()->GetOutputDirectory() == "") || (HeartConfig::Instance()->GetOutputFilenamePrefix() == ""))
242 {
243 EXCEPTION("Either explicitly specify not to print output (call PrintOutput(false)) or specify the output directory and filename prefix");
244 }
245 }
246
247 double end_time = HeartConfig::Instance()->GetSimulationDuration();
248 double pde_time = HeartConfig::Instance()->GetPdeTimeStep();
249
250 /*
251 * MatrixIsConstant stuff requires CONSTANT dt - do some checks to make sure
252 * the TimeStepper won't find non-constant dt.
253 * Note: printing_time does not have to divide end_time, but dt must divide
254 * printing_time and end_time.
255 * HeartConfig checks pde_dt divides printing dt.
256 */
258 if (fabs(end_time - pde_time * round(end_time / pde_time)) > 1e-10)
259 {
260 EXCEPTION("PDE timestep does not seem to divide end time - check parameters");
261 }
262}
263
264template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
266{
267 DistributedVectorFactory* p_factory = mpMesh->GetDistributedVectorFactory();
268 Vec initial_condition = p_factory->CreateVec(PROBLEM_DIM);
269 DistributedVector ic = p_factory->CreateDistributedVector(initial_condition);
270 std::vector<DistributedVector::Stripe> stripe;
271 stripe.reserve(PROBLEM_DIM);
272
273 for (unsigned i = 0; i < PROBLEM_DIM; i++)
274 {
275 stripe.push_back(DistributedVector::Stripe(ic, i));
276 }
277
278 for (DistributedVector::Iterator index = ic.Begin();
279 index != ic.End();
280 ++index)
281 {
282 stripe[0][index] = mpCardiacTissue->GetCardiacCell(index.Global)->GetVoltage();
283 if (PROBLEM_DIM == 2)
284 {
285 stripe[1][index] = 0;
286 }
287 }
288
289 ic.Restore();
290
291 return initial_condition;
292}
293
294template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
296{
297 /*
298 * If this fails the mesh has already been set. We assert rather throw
299 * an exception to avoid a memory leak when checking it throws correctly.
300 */
301 assert(mpMesh == NULL);
302 assert(pMesh != NULL);
303 mAllocatedMemoryForMesh = false;
304 mpMesh = pMesh;
305}
306
307template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
309{
310 mPrintOutput = printOutput;
311}
312
313template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
315{
316 mWriteInfo = writeInfo;
317}
318
319template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
324
325template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
327{
328 return mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(mSolution);
329}
330
331template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
336
337template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
343
344template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
346{
347 if (mpCardiacTissue == NULL)
348 {
349 EXCEPTION("Tissue not yet set up, you may need to call Initialise() before GetTissue().");
350 }
351 return mpCardiacTissue;
352}
353
354template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
356 bool useAdaptivity,
358{
359 if (useAdaptivity)
360 {
361 assert(pController);
362 mpTimeAdaptivityController = pController;
363 }
364 else
365 {
366 mpTimeAdaptivityController = NULL;
367 }
368}
369
370template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
372{
373 PreSolveChecks();
374
375 std::vector<double> additional_stopping_times;
376 SetUpAdditionalStoppingTimes(additional_stopping_times);
377
378 TimeStepper stepper(mCurrentTime,
379 HeartConfig::Instance()->GetSimulationDuration(),
380 HeartConfig::Instance()->GetPrintingTimeStep(),
381 false,
382 additional_stopping_times);
383 // Note that SetUpAdditionalStoppingTimes is a method from the BidomainWithBath class it adds
384 // electrode events into the regular time-stepping
385 // EXCEPTION("Electrode switch on/off events should coincide with printing time steps.");
386
387 if (!mpBoundaryConditionsContainer) // the user didn't supply a bcc
388 {
389 // Set up the default bcc
390 mpDefaultBoundaryConditionsContainer.reset(new BoundaryConditionsContainer<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>);
391 for (unsigned problem_index = 0; problem_index < PROBLEM_DIM; problem_index++)
392 {
393 mpDefaultBoundaryConditionsContainer->DefineZeroNeumannOnMeshBoundary(mpMesh, problem_index);
394 }
395 mpBoundaryConditionsContainer = mpDefaultBoundaryConditionsContainer;
396 }
397
398 assert(mpSolver == NULL);
399 mpSolver = CreateSolver(); // passes mpBoundaryConditionsContainer to solver
400
401 // If we have already run a simulation, use the old solution as initial condition
402 Vec initial_condition;
403 if (mSolution)
404 {
405 initial_condition = mSolution;
406 }
407 else
408 {
409 initial_condition = CreateInitialCondition();
410 }
411
412 std::string progress_reporter_dir;
413
414 if (mPrintOutput)
415 {
416 HeartEventHandler::BeginEvent(HeartEventHandler::WRITE_OUTPUT);
417 bool extending_file = false;
418 try
419 {
420 extending_file = InitialiseWriter();
421 }
422 catch (Exception& e)
423 {
424 delete mpWriter;
425 mpWriter = NULL;
426 delete mpSolver;
427 if (mSolution != initial_condition)
428 {
429 /*
430 * A PETSc Vec is a pointer, so we *don't* need to free the memory if it is
431 * freed somewhere else (e.g. in the destructor). If this is a resumed solution
432 * we set initial_condition = mSolution earlier. mSolution is going to be
433 * cleaned up in the constructor. So, only PetscTools::Destroy( initial_condition ) when
434 * it is not equal to mSolution.
435 */
436 PetscTools::Destroy(initial_condition);
437 }
438 throw e;
439 }
440
441 /*
442 * If we are resuming a simulation (i.e. mSolution already exists) and
443 * we are extending a .h5 file that already exists then there is no need
444 * to write the initial condition to file - it is already there as the
445 * final solution of the previous run.
446 */
447 if (!(mSolution && extending_file))
448 {
449 WriteOneStep(stepper.GetTime(), initial_condition);
450 mpWriter->AdvanceAlongUnlimitedDimension();
451 }
452 HeartEventHandler::EndEvent(HeartEventHandler::WRITE_OUTPUT);
453
454 progress_reporter_dir = HeartConfig::Instance()->GetOutputDirectory();
455 }
456 else
457 {
458 progress_reporter_dir = ""; // progress printed to CHASTE_TEST_OUTPUT
459 }
460 BOOST_FOREACH (boost::shared_ptr<AbstractOutputModifier> p_output_modifier, mOutputModifiers)
462 p_output_modifier->InitialiseAtStart(this->mpMesh->GetDistributedVectorFactory(), this->mpMesh->rGetNodePermutation());
463 p_output_modifier->ProcessSolutionAtTimeStep(stepper.GetTime(), initial_condition, PROBLEM_DIM);
464 }
465
466 /*
467 * Create a progress reporter so users can track how much has gone and
468 * estimate how much time is left. Note this has to be done after the
469 * InitialiseWriter above (if mPrintOutput==true).
470 */
471 ProgressReporter progress_reporter(progress_reporter_dir,
472 mCurrentTime,
473 HeartConfig::Instance()->GetSimulationDuration());
474 progress_reporter.Update(mCurrentTime);
475
476 mpSolver->SetTimeStep(HeartConfig::Instance()->GetPdeTimeStep());
477 if (mpTimeAdaptivityController)
478 {
479 mpSolver->SetTimeAdaptivityController(mpTimeAdaptivityController);
480 }
481
482 while (!stepper.IsTimeAtEnd())
483 {
484 // Solve from now up to the next printing time
485 mpSolver->SetTimes(stepper.GetTime(), stepper.GetNextTime());
486 mpSolver->SetInitialCondition(initial_condition);
487
488 AtBeginningOfTimestep(stepper.GetTime());
489
490 try
491 {
492 try
493 {
494 mSolution = mpSolver->Solve();
496 catch (const Exception& e)
497 {
498#ifndef NDEBUG
500#endif
501 throw e;
502 }
503#ifndef NDEBUG
505#endif
506 }
507 catch (const Exception& e)
508 {
509 // Free memory
510 delete mpSolver;
511 mpSolver = NULL;
512 if (initial_condition != mSolution)
513 {
514 /*
515 * A PETSc Vec is a pointer, so we *don't* need to free the memory if it is
516 * freed somewhere else (e.g. in the destructor). Later, in this while loop
517 * we will set initial_condition = mSolution (or, if this is a resumed solution
518 * it may also have been done when initial_condition was created). mSolution
519 * is going to be cleaned up in the destructor. So, only PetscTools::Destroy()
520 * initial_condition when it is not equal to mSolution (see #1695).
521 */
522 PetscTools::Destroy(initial_condition);
523 }
524
525 // Re-throw
527 CloseFilesAndPostProcess();
528
529 throw e;
530 }
531
532 // Free old initial condition
533 HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
534 PetscTools::Destroy(initial_condition);
535 HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
537 // Initial condition for next loop is current solution
538 initial_condition = mSolution;
539
540 // Update the current time
541 stepper.AdvanceOneTimeStep();
542 mCurrentTime = stepper.GetTime();
543
544 // Print out details at current time if asked for
545 if (mWriteInfo)
546 {
547 HeartEventHandler::BeginEvent(HeartEventHandler::WRITE_OUTPUT);
548 WriteInfo(stepper.GetTime());
549 HeartEventHandler::EndEvent(HeartEventHandler::WRITE_OUTPUT);
550 }
552 BOOST_FOREACH (boost::shared_ptr<AbstractOutputModifier> p_output_modifier, mOutputModifiers)
553 {
554 p_output_modifier->ProcessSolutionAtTimeStep(stepper.GetTime(), mSolution, PROBLEM_DIM);
555 }
556 if (mPrintOutput)
557 {
558 // Writing data out to the file <FilenamePrefix>.dat
559 HeartEventHandler::BeginEvent(HeartEventHandler::WRITE_OUTPUT);
560 WriteOneStep(stepper.GetTime(), mSolution);
561 // Just flags that we've finished a time-step; won't actually 'extend' unless new data is written.
562 mpWriter->AdvanceAlongUnlimitedDimension();
563
564 HeartEventHandler::EndEvent(HeartEventHandler::WRITE_OUTPUT);
566
567 progress_reporter.Update(stepper.GetTime());
568
569 OnEndOfTimestep(stepper.GetTime());
570 }
571
572 // Free solver
573 delete mpSolver;
574 mpSolver = NULL;
575
576 // Close the file that stores voltage values
577 progress_reporter.PrintFinalising();
578 BOOST_FOREACH (boost::shared_ptr<AbstractOutputModifier> p_output_modifier, mOutputModifiers)
579 {
580 p_output_modifier->FinaliseAtEnd();
581 }
582 CloseFilesAndPostProcess();
583 HeartEventHandler::EndEvent(HeartEventHandler::EVERYTHING);
585
586template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
588{
589 // Close files
590 if (!mPrintOutput)
592 // Nothing to do
593 return;
594 }
595 HeartEventHandler::BeginEvent(HeartEventHandler::WRITE_OUTPUT);
596 // If write caching is on, the next line might actually take a significant amount of time.
597 delete mpWriter;
598 mpWriter = NULL;
599 HeartEventHandler::EndEvent(HeartEventHandler::WRITE_OUTPUT);
600
601 FileFinder test_output(HeartConfig::Instance()->GetOutputDirectory(), RelativeTo::ChasteTestOutput);
602
603 /********************************************************************************
604 * Run all post processing.
605 *
606 * The PostProcessingWriter class examines what is requested in HeartConfig and
607 * adds the relevant data to the HDF5 file.
608 * This is converted to different visualizer formats along with the solution
609 * in the DATA_CONVERSION block below.
610 *********************************************************************************/
611
612 HeartEventHandler::BeginEvent(HeartEventHandler::POST_PROC);
613 if (HeartConfig::Instance()->IsPostProcessingRequested())
614 {
616 test_output,
617 HeartConfig::Instance()->GetOutputFilenamePrefix(),
618 "V",
619 mHdf5DataWriterChunkSizeAndAlignment);
620 post_writer.WritePostProcessingFiles();
621 }
622 HeartEventHandler::EndEvent(HeartEventHandler::POST_PROC);
623
624 /********************************************************************************************
625 * Convert HDF5 datasets (solution and postprocessing maps) to different visualizer formats
626 ********************************************************************************************/
627
628 HeartEventHandler::BeginEvent(HeartEventHandler::DATA_CONVERSION);
629 // Only if results files were written and we are outputting all nodes
630 if (mNodesToOutput.empty())
631 {
632 if (HeartConfig::Instance()->GetVisualizeWithMeshalyzer())
633 {
634 // Convert simulation data to Meshalyzer format
636 HeartConfig::Instance()->GetOutputFilenamePrefix(),
637 mpMesh,
638 HeartConfig::Instance()->GetOutputUsingOriginalNodeOrdering(),
639 HeartConfig::Instance()->GetVisualizerOutputPrecision());
640 std::string subdirectory_name = converter.GetSubdirectory();
641 HeartConfig::Instance()->Write(false, subdirectory_name);
642 }
643
644 if (HeartConfig::Instance()->GetVisualizeWithCmgui())
645 {
646 // Convert simulation data to Cmgui format
648 HeartConfig::Instance()->GetOutputFilenamePrefix(),
649 mpMesh,
650 GetHasBath(),
651 HeartConfig::Instance()->GetVisualizerOutputPrecision());
652 std::string subdirectory_name = converter.GetSubdirectory();
653 HeartConfig::Instance()->Write(false, subdirectory_name);
654 }
655
656 if (HeartConfig::Instance()->GetVisualizeWithVtk())
657 {
658 // Convert simulation data to VTK format
659 Hdf5ToVtkConverter<ELEMENT_DIM, SPACE_DIM> converter(test_output,
660 HeartConfig::Instance()->GetOutputFilenamePrefix(),
661 mpMesh,
662 false,
663 HeartConfig::Instance()->GetOutputUsingOriginalNodeOrdering());
664 std::string subdirectory_name = converter.GetSubdirectory();
665 HeartConfig::Instance()->Write(false, subdirectory_name);
666 }
667
668 if (HeartConfig::Instance()->GetVisualizeWithParallelVtk())
669 {
670 // Convert simulation data to parallel VTK (pvtu) format
671 Hdf5ToVtkConverter<ELEMENT_DIM, SPACE_DIM> converter(test_output,
672 HeartConfig::Instance()->GetOutputFilenamePrefix(),
673 mpMesh,
674 true,
675 HeartConfig::Instance()->GetOutputUsingOriginalNodeOrdering());
676 std::string subdirectory_name = converter.GetSubdirectory();
677 HeartConfig::Instance()->Write(false, subdirectory_name);
678 }
679 }
680 HeartEventHandler::EndEvent(HeartEventHandler::DATA_CONVERSION);
681}
682
683template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
685{
686 if (!extending)
687 {
688 if (mNodesToOutput.empty())
689 {
690 //Set writer to output all nodes
691 mpWriter->DefineFixedDimension(mpMesh->GetNumNodes());
692 }
693 else
694 {
695 // Added for #2980
696 if (mpMesh->rGetNodePermutation().size() > 0)
697 {
698 if (HeartConfig::Instance()->GetOutputUsingOriginalNodeOrdering())
699 {
700 EXCEPTION("HeartConfig setting `GetOutputUsingOriginalNodeOrdering` is meaningless when outputting particular nodes in parallel. (Nodes are written with their original indices by default).");
702 std::vector<unsigned> nodes_to_output_permuted(mNodesToOutput.size());
703 for (unsigned i = 0; i < mNodesToOutput.size(); i++)
704 {
705 nodes_to_output_permuted[i] = mpMesh->rGetNodePermutation()[mNodesToOutput[i]];
706 }
707 mpWriter->DefineFixedDimension(mNodesToOutput, nodes_to_output_permuted, mpMesh->GetNumNodes());
708 } else {
709 // Output only the nodes indicated
710 mpWriter->DefineFixedDimension(mNodesToOutput, mNodesToOutput, mpMesh->GetNumNodes());
711 }
713 // mNodeColumnId = mpWriter->DefineVariable("Node", "dimensionless");
714 mVoltageColumnId = mpWriter->DefineVariable("V", "mV");
715
716 // Only used to get an estimate of the # of timesteps below
717 TimeStepper stepper(mCurrentTime,
718 HeartConfig::Instance()->GetSimulationDuration(),
719 HeartConfig::Instance()->GetPrintingTimeStep());
720
721 mpWriter->DefineUnlimitedDimension("Time", "msecs", stepper.EstimateTimeSteps() + 1); // plus one for start and end points
722 }
723 else
724 {
725 mVoltageColumnId = mpWriter->GetVariableByName("V");
726 }
727}
728
729template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
731{
732 mExtraVariablesId.clear();
733 // Check if any extra output variables have been requested
734 if (HeartConfig::Instance()->GetOutputVariablesProvided())
735 {
736 // Get their names in a vector
737 std::vector<std::string> output_variables;
738 HeartConfig::Instance()->GetOutputVariables(output_variables);
739 const unsigned num_vars = output_variables.size();
740 mExtraVariablesId.reserve(num_vars);
741
742 // Loop over them
743 for (unsigned var_index = 0; var_index < num_vars; var_index++)
744 {
745 // Get variable name
746 std::string var_name = output_variables[var_index];
747
748 // Register it (or look it up) in the data writer
749 unsigned column_id;
750 if (extending)
751 {
752 column_id = this->mpWriter->GetVariableByName(var_name);
753 }
754 else
755 {
756 // Difficult to specify the units, as different cell models
757 // at different points in the mesh could be using different units.
758 column_id = this->mpWriter->DefineVariable(var_name, "unknown_units");
759 }
760
761 // Store column id
762 mExtraVariablesId.push_back(column_id);
763 }
764 }
765}
766
767template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
769{
770 // Get the variable names in a vector
771 std::vector<std::string> output_variables;
772 unsigned num_vars = mExtraVariablesId.size();
773 if (num_vars > 0)
774 {
775 HeartConfig::Instance()->GetOutputVariables(output_variables);
776 }
777 assert(output_variables.size() == num_vars);
778
779 // Loop over the requested variables
780 for (unsigned var_index = 0; var_index < num_vars; var_index++)
781 {
782 // Create vector for storing values over the local nodes
783 Vec variable_data = this->mpMesh->GetDistributedVectorFactory()->CreateVec();
784 DistributedVector distributed_var_data = this->mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(variable_data);
785
786 // Loop over the local nodes and gather the data
787 for (DistributedVector::Iterator index = distributed_var_data.Begin();
788 index != distributed_var_data.End();
789 ++index)
790 {
791 // If the region is in the bath
792 if (HeartRegionCode::IsRegionBath(this->mpMesh->GetNode(index.Global)->GetRegion()))
793 {
794 // Then we just pad the output with zeros, user currently needs to find a nice
795 // way to deal with this in processing and visualization.
796 distributed_var_data[index] = 0.0;
797 }
798 else
799 {
800 // Find the variable in the cell model and store its value
801 distributed_var_data[index] = this->mpCardiacTissue->GetCardiacCell(index.Global)->GetAnyVariable(output_variables[var_index], mCurrentTime);
802 }
803 }
804 distributed_var_data.Restore();
805
806 // Write it to disc
807 this->mpWriter->PutVector(mExtraVariablesId[var_index], variable_data);
808
809 PetscTools::Destroy(variable_data);
810 }
811}
812
813template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
815{
816 bool extend_file = (mSolution != NULL);
817
818 // I think this is impossible to trip; certainly it's very difficult!
819 assert(!mpWriter);
820
821 if (extend_file)
822 {
824 + "/" + HeartConfig::Instance()->GetOutputFilenamePrefix() + ".h5",
826 //We are going to test for existence before creating the file.
827 //Therefore we should make sure that this existence test is thread-safe.
828 //(If another process creates the file too early then we may get the wrong answer to the
829 //existence question).
830 PetscTools::Barrier("InitialiseWriter::Extension check");
831 if (!h5_file.Exists())
832 {
833 extend_file = false;
834 }
835 else // if it does exist check that it is sensible to extend it by running from the archive we loaded.
836 {
837 Hdf5DataReader reader(HeartConfig::Instance()->GetOutputDirectory(),
838 HeartConfig::Instance()->GetOutputFilenamePrefix(),
839 true);
840 std::vector<double> times = reader.GetUnlimitedDimensionValues();
841 if (times.back() > mCurrentTime)
842 {
843 EXCEPTION("Attempting to extend " << h5_file.GetAbsolutePath() << " with results from time = " << mCurrentTime << ", but it already contains results up to time = " << times.back() << "."
844 " Calling HeartConfig::Instance()->SetOutputDirectory() before Solve() will direct results elsewhere.");
845 }
846 }
847 PetscTools::Barrier("InitialiseWriter::Extension check");
848 }
849 mpWriter = new Hdf5DataWriter(*mpMesh->GetDistributedVectorFactory(),
852 !extend_file, // don't clear directory if extension requested
853 extend_file,
854 "Data",
855 mUseHdf5DataWriterCache);
856
857 /* If user has specified a chunk size and alignment parameter, pass it
858 * through. We set them to the same value as we think this is the most
859 * likely use case, specifically on striped filesystems where a chunk
860 * should squeeze into a stripe.
861 * Only happens if !extend_file, i.e. we're NOT loading a checkpoint, or
862 * we are loading a checkpoint but the H5 file doesn't exist yet.
863 */
864 if (!extend_file && mHdf5DataWriterChunkSizeAndAlignment)
865 {
866 mpWriter->SetTargetChunkSize(mHdf5DataWriterChunkSizeAndAlignment);
867 mpWriter->SetAlignment(mHdf5DataWriterChunkSizeAndAlignment);
868 }
869
870 // Define columns, or get the variable IDs from the writer
871 DefineWriterColumns(extend_file);
872
873 // Possibility of applying a permutation
874 if (HeartConfig::Instance()->GetOutputUsingOriginalNodeOrdering())
875 {
876 bool success = mpWriter->ApplyPermutation(mpMesh->rGetNodePermutation(), true /*unsafe mode - extending*/);
877 if (success == false)
878 {
879 //It's not really a permutation, so reset
881 }
882 }
883
884 if (!extend_file)
885 {
886 mpWriter->EndDefineMode();
887 }
888
889 return extend_file;
890}
891
892template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
894{
895 mUseHdf5DataWriterCache = useCache;
896}
897
898template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
903
904template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
906{
907 mNodesToOutput = nodesToOutput;
908}
909
910template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
912{
913 if ((HeartConfig::Instance()->GetOutputDirectory() == "") || (HeartConfig::Instance()->GetOutputFilenamePrefix() == ""))
914 {
915 EXCEPTION("Data reader invalid as data writer cannot be initialised");
916 }
917 return Hdf5DataReader(HeartConfig::Instance()->GetOutputDirectory(), HeartConfig::Instance()->GetOutputFilenamePrefix());
918}
919
920template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
925
926template <unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
930
931// Explicit instantiation
932
933// Monodomain
939
940// Bidomain
944
945// Extended Bidomain
#define EXCEPTION(message)
#define NEVER_REACHED
void SetWriteInfo(bool writeInfo=true)
void SetOutputNodes(std::vector< unsigned > &rNodesToOutput)
void DefineExtraVariablesWriterColumns(bool extending)
DistributedVector GetSolutionDistributedVector()
void SetHdf5DataWriterTargetChunkSizeAndAlignment(hsize_t size)
void SetUseTimeAdaptivityController(bool useAdaptivity, AbstractTimeAdaptivityController *pController=NULL)
void SetMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh)
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > & rGetMesh()
void PrintOutput(bool rPrintOutput)
void SetBoundaryConditionsContainer(BccType pBcc)
AbstractCardiacCellFactory< ELEMENT_DIM, SPACE_DIM > * mpCellFactory
std::vector< unsigned > mNodesToOutput
AbstractCardiacTissue< ELEMENT_DIM, SPACE_DIM > * GetTissue()
void SetUseHdf5DataWriterCache(bool useCache=true)
virtual void DefineWriterColumns(bool extending)
DistributedVector CreateDistributedVector(Vec vec, bool readOnly=false)
std::string GetAbsolutePath() const
bool Exists() const
std::vector< double > GetUnlimitedDimensionValues()
void GetSheetDimensions(c_vector< double, 2 > &sheetDimensions) const
double GetPdeTimeStep() const
void GetOutputVariables(std::vector< std::string > &rOutputVariables) const
double GetSimulationDuration() const
void SetOutputUsingOriginalNodeOrdering(bool useOriginal)
std::string GetOutputFilenamePrefix() const
void GetSlabDimensions(c_vector< double, 3 > &slabDimensions) const
void Write(bool useArchiveLocationInfo=false, std::string subfolderName="output")
double GetInterNodeSpace() const
std::string GetOutputDirectory() const
void GetFibreLength(c_vector< double, 1 > &fibreLength) const
static HeartConfig * Instance()
static bool IsRegionBath(HeartRegionType regionId)
static std::string GetChasteTestOutputDirectory()
static void Destroy(Vec &rVec)
static void Barrier(const std::string callerId="")
static bool IsParallel()
static void ReplicateException(bool flag)
bool IsTimeAtEnd() const
double GetTime() const
void AdvanceOneTimeStep()
double GetNextTime() const
unsigned EstimateTimeSteps() const