Chaste  Release::2018.1
OperatorSplittingMonodomainSolver.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, 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 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 #include "OperatorSplittingMonodomainSolver.hpp"
37 
38 
39 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
41 {
42  assert(this->mpLinearSystem->rGetLhsMatrix() != NULL);
43  assert(this->mpLinearSystem->rGetRhsVector() != NULL);
44 
46  // set up LHS matrix (and mass matrix)
48  if (computeMatrix)
49  {
50  mpMonodomainAssembler->SetMatrixToAssemble(this->mpLinearSystem->rGetLhsMatrix());
51  mpMonodomainAssembler->AssembleMatrix();
52 
53  MassMatrixAssembler<ELEMENT_DIM,SPACE_DIM> mass_matrix_assembler(this->mpMesh, HeartConfig::Instance()->GetUseMassLumping());
54  mass_matrix_assembler.SetMatrixToAssemble(mMassMatrix);
55  mass_matrix_assembler.Assemble();
56 
57  this->mpLinearSystem->FinaliseLhsMatrix();
58  PetscMatTools::Finalise(mMassMatrix);
59  }
60 
61  HeartEventHandler::BeginEvent(HeartEventHandler::ASSEMBLE_RHS);
62 
64  // Set up z in b=Mz
66  DistributedVectorFactory* p_factory = this->mpMesh->GetDistributedVectorFactory();
67  // dist stripe for the current Voltage
68  DistributedVector distributed_current_solution = p_factory->CreateDistributedVector(currentSolution);
69  // dist stripe for z (return value)
70  DistributedVector dist_vec_matrix_based = p_factory->CreateDistributedVector(mVecForConstructingRhs);
71 
73  double Cm = HeartConfig::Instance()->GetCapacitance();
74 
75  for (DistributedVector::Iterator index = dist_vec_matrix_based.Begin();
76  index!= dist_vec_matrix_based.End();
77  ++index)
78  {
79  double V = distributed_current_solution[index];
80  // in the main solver, the nodal ionic current and stimuli is computed and used.
81  // However in operator splitting, this part of the solve is diffusion only, no reaction terms
82  //double F = - Am*this->mpMonodomainTissue->rGetIionicCacheReplicated()[index.Global]
83  // - this->mpMonodomainTissue->rGetIntracellularStimulusCacheReplicated()[index.Global];
84 
85  dist_vec_matrix_based[index] = Am*Cm*V*PdeSimulationTime::GetPdeTimeStepInverse();
86  }
87  dist_vec_matrix_based.Restore();
88 
90  // b = Mz
92  MatMult(mMassMatrix, mVecForConstructingRhs, this->mpLinearSystem->rGetRhsVector());
93 
94  // assembling RHS is not finished yet, as Neumann bcs are added below, but
95  // the event will be begun again inside mpMonodomainAssembler->AssembleVector();
96  HeartEventHandler::EndEvent(HeartEventHandler::ASSEMBLE_RHS);
97 
99  // apply Neumann boundary conditions
101  mpNeumannSurfaceTermsAssembler->SetVectorToAssemble(this->mpLinearSystem->rGetRhsVector(), false/*don't zero vector!*/);
102  mpNeumannSurfaceTermsAssembler->AssembleVector();
103 
104  // finalise
105  this->mpLinearSystem->FinaliseRhsVector();
106 }
107 
108 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
110 {
111  double time = PdeSimulationTime::GetTime();
112  double dt = PdeSimulationTime::GetPdeTimeStep();
113  mpMonodomainTissue->SolveCellSystems(currentSolution, time, time+dt/2.0, true);
114 }
115 
116 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
118 {
119  // solve cell models for second half timestep
120  double time = PdeSimulationTime::GetTime();
121  double dt = PdeSimulationTime::GetPdeTimeStep();
122  mpMonodomainTissue->SolveCellSystems(currentSolution, time + dt/2, PdeSimulationTime::GetNextTime(), true);
123 }
124 
125 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
127 {
128  if (this->mpLinearSystem != NULL)
129  {
130  return;
131  }
132 
133  // call base class version...
135 
136  //..then do a bit extra
137  if (HeartConfig::Instance()->GetUseAbsoluteTolerance())
138  {
139  this->mpLinearSystem->SetAbsoluteTolerance(HeartConfig::Instance()->GetAbsoluteTolerance());
140  }
141  else
142  {
144  // re-implement when needed
145  //this->mpLinearSystem->SetRelativeTolerance(HeartConfig::Instance()->GetRelativeTolerance());
146  }
147 
148  this->mpLinearSystem->SetKspType(HeartConfig::Instance()->GetKSPSolver());
149  this->mpLinearSystem->SetPcType(HeartConfig::Instance()->GetKSPPreconditioner());
150  this->mpLinearSystem->SetMatrixIsSymmetric(true);
151  this->mpLinearSystem->SetUseFixedNumberIterations(HeartConfig::Instance()->GetUseFixedNumberIterationsLinearSolver(), HeartConfig::Instance()->GetEvaluateNumItsEveryNSolves());
152 
153  // initialise matrix-based RHS vector and matrix, and use the linear
154  // system rhs as a template
155  Vec& r_template = this->mpLinearSystem->rGetRhsVector();
156  VecDuplicate(r_template, &mVecForConstructingRhs);
157  PetscInt ownership_range_lo;
158  PetscInt ownership_range_hi;
159  VecGetOwnershipRange(r_template, &ownership_range_lo, &ownership_range_hi);
160  PetscInt local_size = ownership_range_hi - ownership_range_lo;
161  PetscTools::SetupMat(mMassMatrix, this->mpMesh->GetNumNodes(), this->mpMesh->GetNumNodes(),
162  this->mpMesh->CalculateMaximumNodeConnectivityPerProcess(),
163  local_size, local_size);
164 }
165 
166 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
171  : AbstractDynamicLinearPdeSolver<ELEMENT_DIM,SPACE_DIM,1>(pMesh),
172  mpBoundaryConditions(pBoundaryConditions),
173  mpMonodomainTissue(pTissue)
174 {
175  assert(pTissue);
176  assert(pBoundaryConditions);
177  this->mMatrixIsConstant = true;
178 
181 
182  // Tell tissue there's no need to replicate ionic caches
183  pTissue->SetCacheReplication(false);
184  mVecForConstructingRhs = NULL;
185 }
186 
187 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
189 {
190  delete mpMonodomainAssembler;
191  delete mpNeumannSurfaceTermsAssembler;
192 
193  if (mVecForConstructingRhs)
194  {
195  PetscTools::Destroy(mVecForConstructingRhs);
196  PetscTools::Destroy(mMassMatrix);
197  }
198 }
199 
200 // Explicit instantiation
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
DistributedVector CreateDistributedVector(Vec vec, bool readOnly=false)
virtual void InitialiseForSolve(Vec initialSolution=nullptr)
NaturalNeumannSurfaceTermAssembler< ELEMENT_DIM, SPACE_DIM, 1 > * mpNeumannSurfaceTermsAssembler
#define NEVER_REACHED
Definition: Exception.hpp:206
OperatorSplittingMonodomainSolver(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, MonodomainTissue< ELEMENT_DIM, SPACE_DIM > *pTissue, BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, 1 > *pBoundaryConditions)
MonodomainTissue< ELEMENT_DIM, SPACE_DIM > * mpMonodomainTissue
void SetupLinearSystem(Vec currentSolution, bool computeMatrix)
static double GetNextTime()
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)
Definition: PetscTools.cpp:268
static double GetPdeTimeStepInverse()
static void Destroy(Vec &rVec)
Definition: PetscTools.hpp:352
static double GetPdeTimeStep()
void SetCacheReplication(bool doCacheReplication)
static void Finalise(Mat matrix)
double GetCapacitance() const
double GetSurfaceAreaToVolumeRatio() const
void SetMatrixToAssemble(Mat &rMatToAssemble, bool zeroMatrixBeforeAssembly=true)
static HeartConfig * Instance()
static double GetTime()
MonodomainAssembler< ELEMENT_DIM, SPACE_DIM > * mpMonodomainAssembler