#include <OperatorSplittingMonodomainSolver.hpp>
Inherits AbstractDynamicLinearPdeSolver< ELEMENT_DIM, SPACE_DIM, 1 >.
Public Member Functions | |
void | InitialiseForSolve (Vec initialSolution) |
OperatorSplittingMonodomainSolver (AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, MonodomainTissue< ELEMENT_DIM, SPACE_DIM > *pTissue, BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, 1 > *pBoundaryConditions, unsigned numQuadPoints=2) | |
~OperatorSplittingMonodomainSolver () | |
Private Member Functions | |
void | SetupLinearSystem (Vec currentSolution, bool computeMatrix) |
void | PrepareForSetupLinearSystem (Vec currentSolution) |
void | FollowingSolveLinearSystem (Vec currentSolution) |
Private Attributes | |
BoundaryConditionsContainer < ELEMENT_DIM, SPACE_DIM, 1 > * | mpBoundaryConditions |
MonodomainTissue< ELEMENT_DIM, SPACE_DIM > * | mpMonodomainTissue |
MonodomainAssembler < ELEMENT_DIM, SPACE_DIM > * | mpMonodomainAssembler |
NaturalNeumannSurfaceTermAssembler < ELEMENT_DIM, SPACE_DIM, 1 > * | mpNeumannSurfaceTermsAssembler |
unsigned | mNumQuadPoints |
Mat | mMassMatrix |
Vec | mVecForConstructingRhs |
A monodomain solver that uses Strang operator splitting of the diffusion (conductivity) term and the reaction (ionic current) term, instead of solving the full reaction-diffusion PDE. This does NOT refer to operator splitting of the two PDEs in the bidomain equations. For details see for example Sundnes et al "Computing the Electrical Activity of the Heart".
The algorithm is, for solving from t=T to T+dt.
(i) Solve ODEs dV/dt = Iionic for t=T to T+dt/2 [giving updated V (internally, and in solution vector) and updated state variables] (ii) Solve PDE dV/dt = div sigma grad V for t=T to dt [using V from step i, --> updated V] (iii) Solve ODEs dV/dt = Iionic for t=T+dt/2 to T+dt [using V from step ii, --> final V]
Notes (a) Stages (iii) and (i) can normally be solved together in one go, except just before/after printing the voltage to file. However for simplicity of code this has not been implemented (b) Therefore, the effective ODE timestep will be: min(ode_dt, pde_dt/2), where ode_dt and pde_dt are those given via HeartConfig. (c) This solver is FOR COMPARING ACCURACY, NOT PERFORMANCE. It has not been optimised and may or may not perform well in parallel. (d) We don't implement the simpler form of operator splitting, Godunov splitting, where the ODEs are solved for one timestep and the PDEs are solved for one timestep, since this is formally equivalent to the default implementation where the ionic current is interpolated from the nodal values (ie ICI - see ICI/SVI discussion in documentation)
Definition at line 65 of file OperatorSplittingMonodomainSolver.hpp.
OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::OperatorSplittingMonodomainSolver | ( | AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * | pMesh, | |
MonodomainTissue< ELEMENT_DIM, SPACE_DIM > * | pTissue, | |||
BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, 1 > * | pBoundaryConditions, | |||
unsigned | numQuadPoints = 2 | |||
) | [inline] |
Constructor
pMesh | pointer to the mesh | |
pTissue | pointer to the tissue | |
pBoundaryConditions | pointer to the boundary conditions | |
numQuadPoints | number of quadrature points (defaults to 2) |
Definition at line 164 of file OperatorSplittingMonodomainSolver.cpp.
References AbstractDynamicLinearPdeSolver< ELEMENT_DIM, SPACE_DIM, 1 >::mMatrixIsConstant, OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mNumQuadPoints, AbstractLinearPdeSolver< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM >::mpMesh, OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mpMonodomainAssembler, OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mpMonodomainTissue, OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mpNeumannSurfaceTermsAssembler, OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mVecForConstructingRhs, and AbstractCardiacTissue< ELEMENT_DIM, SPACE_DIM >::SetCacheReplication().
OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::~OperatorSplittingMonodomainSolver | ( | ) | [inline] |
Destructor
Definition at line 187 of file OperatorSplittingMonodomainSolver.cpp.
References OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mMassMatrix, OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mpMonodomainAssembler, OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mpNeumannSurfaceTermsAssembler, and OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mVecForConstructingRhs.
void OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::FollowingSolveLinearSystem | ( | Vec | currentSolution | ) | [inline, private, virtual] |
Called after solving the linear system, used to solve the cell models for second half timestep (step (iii) above)
currentSolution | the latest solution vector (ie the solution of the linear system). |
Reimplemented from AbstractLinearPdeSolver< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM >.
Definition at line 111 of file OperatorSplittingMonodomainSolver.cpp.
References PdeSimulationTime::GetPdeTimeStep(), PdeSimulationTime::GetTime(), and OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mpMonodomainTissue.
void OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::InitialiseForSolve | ( | Vec | initialSolution | ) | [inline, virtual] |
Overloaded InitialiseForSolve() which calls base version but also initialises mMassMatrix and mVecForConstructingRhs.
initialSolution | initial solution |
Reimplemented from AbstractLinearPdeSolver< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM >.
Definition at line 121 of file OperatorSplittingMonodomainSolver.cpp.
References HeartConfig::Instance(), OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mMassMatrix, AbstractLinearPdeSolver< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM >::mpLinearSystem, AbstractLinearPdeSolver< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM >::mpMesh, OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mVecForConstructingRhs, NEVER_REACHED, LinearSystem::rGetRhsVector(), LinearSystem::SetAbsoluteTolerance(), LinearSystem::SetKspType(), LinearSystem::SetMatrixIsSymmetric(), LinearSystem::SetPcType(), PetscTools::SetupMat(), and LinearSystem::SetUseFixedNumberIterations().
void OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::PrepareForSetupLinearSystem | ( | Vec | currentSolution | ) | [inline, private, virtual] |
Called before setting up the linear system, used to solve the cell models for first half timestep (step (i) above)
currentSolution | the latest solution vector |
Reimplemented from AbstractLinearPdeSolver< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM >.
Definition at line 102 of file OperatorSplittingMonodomainSolver.cpp.
References PdeSimulationTime::GetPdeTimeStep(), PdeSimulationTime::GetTime(), and OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mpMonodomainTissue.
void OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::SetupLinearSystem | ( | Vec | currentSolution, | |
bool | computeMatrix | |||
) | [inline, private, virtual] |
Implementation of SetupLinearSystem() which uses the assembler to compute the LHS matrix, but sets up the RHS vector using the mass-matrix (constructed using a separate assembler) multiplied by a vector
currentSolution | Solution at current time | |
computeMatrix | Whether to compute the matrix of the linear system |
Implements AbstractLinearPdeSolver< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM >.
Definition at line 33 of file OperatorSplittingMonodomainSolver.cpp.
References AbstractFeAssemblerInterface< CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX >::Assemble(), AbstractFeAssemblerInterface< CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX >::AssembleVector(), DistributedVector::Begin(), GenericEventHandler< 16, HeartEventHandler >::BeginEvent(), DistributedVectorFactory::CreateDistributedVector(), DistributedVector::End(), GenericEventHandler< 16, HeartEventHandler >::EndEvent(), PetscMatTools::Finalise(), LinearSystem::FinaliseLhsMatrix(), LinearSystem::FinaliseRhsVector(), HeartConfig::GetCapacitance(), PdeSimulationTime::GetPdeTimeStepInverse(), HeartConfig::GetSurfaceAreaToVolumeRatio(), HeartConfig::Instance(), OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mMassMatrix, AbstractLinearPdeSolver< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM >::mpLinearSystem, AbstractLinearPdeSolver< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM >::mpMesh, OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mpMonodomainAssembler, OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mpNeumannSurfaceTermsAssembler, OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mVecForConstructingRhs, DistributedVector::Restore(), LinearSystem::rGetLhsMatrix(), LinearSystem::rGetRhsVector(), AbstractFeAssemblerInterface< CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX >::SetMatrixToAssemble(), and AbstractFeAssemblerInterface< CAN_ASSEMBLE_VECTOR, CAN_ASSEMBLE_MATRIX >::SetVectorToAssemble().
Mat OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mMassMatrix [private] |
The mass matrix, used to computing the RHS vector
Definition at line 88 of file OperatorSplittingMonodomainSolver.hpp.
Referenced by OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::InitialiseForSolve(), OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::SetupLinearSystem(), and OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::~OperatorSplittingMonodomainSolver().
unsigned OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mNumQuadPoints [private] |
Number of quadrature points per dimension (only saved so it can be passed to the assembler
Definition at line 85 of file OperatorSplittingMonodomainSolver.hpp.
Referenced by OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::OperatorSplittingMonodomainSolver().
BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,1>* OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mpBoundaryConditions [private] |
Boundary conditions
Definition at line 70 of file OperatorSplittingMonodomainSolver.hpp.
MonodomainAssembler<ELEMENT_DIM,SPACE_DIM>* OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mpMonodomainAssembler [private] |
The monodomain assembler, used to set up the LHS matrix
Definition at line 76 of file OperatorSplittingMonodomainSolver.hpp.
Referenced by OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::OperatorSplittingMonodomainSolver(), OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::SetupLinearSystem(), and OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::~OperatorSplittingMonodomainSolver().
MonodomainTissue<ELEMENT_DIM,SPACE_DIM>* OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mpMonodomainTissue [private] |
Monodomain tissue class (collection of cells, and conductivities)
Definition at line 73 of file OperatorSplittingMonodomainSolver.hpp.
Referenced by OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::FollowingSolveLinearSystem(), OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::OperatorSplittingMonodomainSolver(), and OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::PrepareForSetupLinearSystem().
NaturalNeumannSurfaceTermAssembler<ELEMENT_DIM,SPACE_DIM,1>* OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mpNeumannSurfaceTermsAssembler [private] |
Assembler for surface integrals coming from any non-zero Neumann boundary conditions
Definition at line 79 of file OperatorSplittingMonodomainSolver.hpp.
Referenced by OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::OperatorSplittingMonodomainSolver(), OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::SetupLinearSystem(), and OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::~OperatorSplittingMonodomainSolver().
Vec OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::mVecForConstructingRhs [private] |
The vector multiplied by the mass matrix. Ie, if the linear system to be solved is Ax=b, this vector is z where b=Mz.
In the normal solver this is equal to alpha*V + beta*Iionic + gamma*Istim, here there is no Iionic term.
Definition at line 97 of file OperatorSplittingMonodomainSolver.hpp.
Referenced by OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::InitialiseForSolve(), OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::OperatorSplittingMonodomainSolver(), OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::SetupLinearSystem(), and OperatorSplittingMonodomainSolver< ELEMENT_DIM, SPACE_DIM >::~OperatorSplittingMonodomainSolver().