Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
MonodomainSolver.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 "MonodomainSolver.hpp"
37#include "MassMatrixAssembler.hpp"
38#include "PetscMatTools.hpp"
39
40
41template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
42void MonodomainSolver<ELEMENT_DIM,SPACE_DIM>::SetupLinearSystem(Vec currentSolution, bool computeMatrix)
43{
44 assert(this->mpLinearSystem->rGetLhsMatrix() != NULL);
45 assert(this->mpLinearSystem->rGetRhsVector() != NULL);
46
47
49 // set up LHS matrix (and mass matrix)
51 if (computeMatrix)
52 {
53 mpMonodomainAssembler->SetMatrixToAssemble(this->mpLinearSystem->rGetLhsMatrix());
54 mpMonodomainAssembler->AssembleMatrix();
55
56 MassMatrixAssembler<ELEMENT_DIM,SPACE_DIM> mass_matrix_assembler(this->mpMesh, HeartConfig::Instance()->GetUseMassLumping());
57 mass_matrix_assembler.SetMatrixToAssemble(mMassMatrix);
58 mass_matrix_assembler.Assemble();
59
60 this->mpLinearSystem->FinaliseLhsMatrix();
61 PetscMatTools::Finalise(mMassMatrix);
62
63 if (HeartConfig::Instance()->GetUseMassLumpingForPrecond() && !HeartConfig::Instance()->GetUseMassLumping())
64 {
65 this->mpLinearSystem->SetPrecondMatrixIsDifferentFromLhs();
66
67 MonodomainAssembler<ELEMENT_DIM,SPACE_DIM> lumped_mass_assembler(this->mpMesh,this->mpMonodomainTissue);
68 lumped_mass_assembler.SetMatrixToAssemble(this->mpLinearSystem->rGetPrecondMatrix());
69
71 lumped_mass_assembler.AssembleMatrix();
73
74 this->mpLinearSystem->FinalisePrecondMatrix();
75 }
76 }
77
78 HeartEventHandler::BeginEvent(HeartEventHandler::ASSEMBLE_RHS);
79
81 // Set up z in b=Mz
83 DistributedVectorFactory* p_factory = this->mpMesh->GetDistributedVectorFactory();
84 // dist stripe for the current Voltage
85 DistributedVector distributed_current_solution = p_factory->CreateDistributedVector(currentSolution);
86 // dist stripe for z (return value)
87 DistributedVector dist_vec_matrix_based = p_factory->CreateDistributedVector(mVecForConstructingRhs);
88
91
92 for (DistributedVector::Iterator index = dist_vec_matrix_based.Begin();
93 index!= dist_vec_matrix_based.End();
94 ++index)
95 {
96 double V = distributed_current_solution[index];
97 double F = - Am*this->mpMonodomainTissue->rGetIionicCacheReplicated()[index.Global]
98 - this->mpMonodomainTissue->rGetIntracellularStimulusCacheReplicated()[index.Global];
99
100 dist_vec_matrix_based[index] = Am*Cm*V*PdeSimulationTime::GetPdeTimeStepInverse() + F;
101 }
102 dist_vec_matrix_based.Restore();
103
105 // b = Mz
107 MatMult(mMassMatrix, mVecForConstructingRhs, this->mpLinearSystem->rGetRhsVector());
108
109 // assembling RHS is not finished yet, as Neumann bcs are added below, but
110 // the event will be begun again inside mpMonodomainAssembler->AssembleVector();
111 HeartEventHandler::EndEvent(HeartEventHandler::ASSEMBLE_RHS);
112
114 // apply Neumann boundary conditions
116 mpNeumannSurfaceTermsAssembler->SetVectorToAssemble(this->mpLinearSystem->rGetRhsVector(), false/*don't zero vector!*/);
117 mpNeumannSurfaceTermsAssembler->AssembleVector();
118
120 // apply correction term
122 if (mpMonodomainCorrectionTermAssembler)
123 {
124 mpMonodomainCorrectionTermAssembler->SetVectorToAssemble(this->mpLinearSystem->rGetRhsVector(), false/*don't zero vector!*/);
125 // don't need to set current solution
126 mpMonodomainCorrectionTermAssembler->AssembleVector();
127 }
128
129 // finalise
130 this->mpLinearSystem->FinaliseRhsVector();
131}
132
133template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
135{
136 if (this->mpLinearSystem != NULL)
137 {
138 return;
139 }
140
141 // call base class version...
143
144 //..then do a bit extra
145 if (HeartConfig::Instance()->GetUseAbsoluteTolerance())
146 {
147 this->mpLinearSystem->SetAbsoluteTolerance(HeartConfig::Instance()->GetAbsoluteTolerance());
148 }
149 else
150 {
151 this->mpLinearSystem->SetRelativeTolerance(HeartConfig::Instance()->GetRelativeTolerance());
152 }
153
154 this->mpLinearSystem->SetKspType(HeartConfig::Instance()->GetKSPSolver());
155 this->mpLinearSystem->SetPcType(HeartConfig::Instance()->GetKSPPreconditioner());
156 this->mpLinearSystem->SetMatrixIsSymmetric(true);
157 this->mpLinearSystem->SetUseFixedNumberIterations(HeartConfig::Instance()->GetUseFixedNumberIterationsLinearSolver(), HeartConfig::Instance()->GetEvaluateNumItsEveryNSolves());
158
159 // initialise matrix-based RHS vector and matrix, and use the linear
160 // system rhs as a template
161 Vec& r_template = this->mpLinearSystem->rGetRhsVector();
162 VecDuplicate(r_template, &mVecForConstructingRhs);
163 PetscInt ownership_range_lo;
164 PetscInt ownership_range_hi;
165 VecGetOwnershipRange(r_template, &ownership_range_lo, &ownership_range_hi);
166 PetscInt local_size = ownership_range_hi - ownership_range_lo;
167 PetscTools::SetupMat(mMassMatrix, this->mpMesh->GetNumNodes(), this->mpMesh->GetNumNodes(),
168 this->mpMesh->CalculateMaximumNodeConnectivityPerProcess(),
169 local_size, local_size);
170}
171
172template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
174{
175 // solve cell models
176 mpMonodomainTissue->SolveCellSystems(currentSolution, PdeSimulationTime::GetTime(), PdeSimulationTime::GetNextTime());
177}
178
179template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
184 : AbstractDynamicLinearPdeSolver<ELEMENT_DIM,SPACE_DIM,1>(pMesh),
185 mpMonodomainTissue(pTissue),
186 mpBoundaryConditions(pBoundaryConditions)
187{
188 assert(pTissue);
189 assert(pBoundaryConditions);
190 this->mMatrixIsConstant = true;
191
194
195
196 // Tell tissue there's no need to replicate ionic caches
197 pTissue->SetCacheReplication(false);
199
200 if (HeartConfig::Instance()->GetUseStateVariableInterpolation())
201 {
204 //We are going to need those caches after all
205 pTissue->SetCacheReplication(true);
206 }
207 else
208 {
210 }
211}
212
213template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
215{
216 delete mpMonodomainAssembler;
217 delete mpNeumannSurfaceTermsAssembler;
218
219 if (mVecForConstructingRhs)
220 {
221 PetscTools::Destroy(mVecForConstructingRhs);
222 PetscTools::Destroy(mMassMatrix);
223 }
224
225 if (mpMonodomainCorrectionTermAssembler)
226 {
227 delete mpMonodomainCorrectionTermAssembler;
228 }
229}
230
231// Explicit instantiation
232template class MonodomainSolver<1,1>;
233template class MonodomainSolver<1,2>;
234template class MonodomainSolver<1,3>;
235template class MonodomainSolver<2,2>;
236template class MonodomainSolver<3,3>;
void SetCacheReplication(bool doCacheReplication)
void SetMatrixToAssemble(Mat &rMatToAssemble, bool zeroMatrixBeforeAssembly=true)
virtual void InitialiseForSolve(Vec initialSolution=nullptr)
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
DistributedVector CreateDistributedVector(Vec vec, bool readOnly=false)
double GetCapacitance() const
double GetSurfaceAreaToVolumeRatio() const
void SetUseMassLumping(bool useMassLumping=true)
static HeartConfig * Instance()
virtual void InitialiseForSolve(Vec initialSolution)
MonodomainAssembler< ELEMENT_DIM, SPACE_DIM > * mpMonodomainAssembler
void PrepareForSetupLinearSystem(Vec currentSolution)
MonodomainCorrectionTermAssembler< ELEMENT_DIM, SPACE_DIM > * mpMonodomainCorrectionTermAssembler
MonodomainSolver(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, MonodomainTissue< ELEMENT_DIM, SPACE_DIM > *pTissue, BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, 1 > *pBoundaryConditions)
void SetupLinearSystem(Vec currentSolution, bool computeMatrix)
NaturalNeumannSurfaceTermAssembler< ELEMENT_DIM, SPACE_DIM, 1 > * mpNeumannSurfaceTermsAssembler
MonodomainTissue< ELEMENT_DIM, SPACE_DIM > * mpMonodomainTissue
static double GetPdeTimeStepInverse()
static double GetTime()
static double GetNextTime()
static void Finalise(Mat matrix)
static void Destroy(Vec &rVec)
static void SetupMat(Mat &rMat, int numRows, int numColumns, unsigned rowPreallocation, int numLocalRows=PETSC_DECIDE, int numLocalColumns=PETSC_DECIDE, bool ignoreOffProcEntries=true, bool newAllocationError=true)