Chaste  Release::2018.1
BidomainSolver.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 
37 #include "BidomainSolver.hpp"
38 #include "BidomainAssembler.hpp"
39 #include "BidomainWithBathAssembler.hpp"
40 #include "PetscMatTools.hpp"
41 
42 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
44 {
45  if (this->mpLinearSystem != NULL)
46  {
47  return;
48  }
50 
51  // initialise matrix-based RHS vector and matrix, and use the linear
52  // system rhs as a template
53  Vec& r_template = this->mpLinearSystem->rGetRhsVector();
54  VecDuplicate(r_template, &mVecForConstructingRhs);
55  PetscInt ownership_range_lo;
56  PetscInt ownership_range_hi;
57  VecGetOwnershipRange(r_template, &ownership_range_lo, &ownership_range_hi);
58  PetscInt local_size = ownership_range_hi - ownership_range_lo;
59  PetscTools::SetupMat(mMassMatrix, 2*this->mpMesh->GetNumNodes(), 2*this->mpMesh->GetNumNodes(),
60  2*this->mpMesh->CalculateMaximumNodeConnectivityPerProcess(),
61  local_size, local_size);
62 }
63 
64 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
66  Vec currentSolution,
67  bool computeMatrix)
68 {
69  assert(this->mpLinearSystem->rGetLhsMatrix() != NULL);
70  assert(this->mpLinearSystem->rGetRhsVector() != NULL);
71  assert(currentSolution != NULL);
72 
73 
75  // set up LHS matrix (and mass matrix)
77  if (computeMatrix)
78  {
79  mpBidomainAssembler->SetMatrixToAssemble(this->mpLinearSystem->rGetLhsMatrix());
80  mpBidomainAssembler->AssembleMatrix();
81 
82  // the BidomainMassMatrixAssembler deals with the mass matrix
83  // for both bath and nonbath problems
84  assert(SPACE_DIM==ELEMENT_DIM);
85  BidomainMassMatrixAssembler<SPACE_DIM> mass_matrix_assembler(this->mpMesh);
86  mass_matrix_assembler.SetMatrixToAssemble(mMassMatrix);
87  mass_matrix_assembler.Assemble();
88 
89  this->mpLinearSystem->SwitchWriteModeLhsMatrix();
90  PetscMatTools::Finalise(mMassMatrix);
91  }
92 
93 
94  HeartEventHandler::BeginEvent(HeartEventHandler::ASSEMBLE_RHS);
95 
97  // Set up z in b=Mz
99  DistributedVectorFactory* p_factory = this->mpMesh->GetDistributedVectorFactory();
100 
101  // dist stripe for the current Voltage
102  DistributedVector distributed_current_solution = p_factory->CreateDistributedVector(currentSolution);
103  DistributedVector::Stripe distributed_current_solution_vm(distributed_current_solution, 0);
104 
105  // dist stripe for z
106  DistributedVector dist_vec_matrix_based = p_factory->CreateDistributedVector(mVecForConstructingRhs);
107  DistributedVector::Stripe dist_vec_matrix_based_vm(dist_vec_matrix_based, 0);
108  DistributedVector::Stripe dist_vec_matrix_based_phie(dist_vec_matrix_based, 1);
109 
111  double Cm = HeartConfig::Instance()->GetCapacitance();
112 
113  if (!(this->mBathSimulation))
114  {
115  for (DistributedVector::Iterator index = dist_vec_matrix_based.Begin();
116  index!= dist_vec_matrix_based.End();
117  ++index)
118  {
119  double V = distributed_current_solution_vm[index];
120  double F = - Am*this->mpBidomainTissue->rGetIionicCacheReplicated()[index.Global]
121  - this->mpBidomainTissue->rGetIntracellularStimulusCacheReplicated()[index.Global];
122 
123  dist_vec_matrix_based_vm[index] = Am*Cm*V*PdeSimulationTime::GetPdeTimeStepInverse() + F;
124  dist_vec_matrix_based_phie[index] = 0.0;
125  }
126  }
127  else
128  {
129  for (DistributedVector::Iterator index = dist_vec_matrix_based.Begin();
130  index!= dist_vec_matrix_based.End();
131  ++index)
132  {
133 
134  if (!HeartRegionCode::IsRegionBath( this->mpMesh->GetNode(index.Global)->GetRegion()))
135  {
136  double V = distributed_current_solution_vm[index];
137  double F = - Am*this->mpBidomainTissue->rGetIionicCacheReplicated()[index.Global]
138  - this->mpBidomainTissue->rGetIntracellularStimulusCacheReplicated()[index.Global];
139 
140  dist_vec_matrix_based_vm[index] = Am*Cm*V*PdeSimulationTime::GetPdeTimeStepInverse() + F;
141  }
142  else
143  {
144  dist_vec_matrix_based_vm[index] = 0.0;
145  }
146 
147  dist_vec_matrix_based_phie[index] = 0.0;
148 
149  }
150  }
151 
152  dist_vec_matrix_based.Restore();
153 
155  // b = Mz
157  MatMult(mMassMatrix, mVecForConstructingRhs, this->mpLinearSystem->rGetRhsVector());
158 
159  // assembling RHS is not finished yet, as Neumann bcs are added below, but
160  // the event will be begun again inside mpBidomainAssembler->AssembleVector();
161  HeartEventHandler::EndEvent(HeartEventHandler::ASSEMBLE_RHS);
162 
163 
165  // apply Neumann boundary conditions
167  mpBidomainNeumannSurfaceTermAssembler->ResetBoundaryConditionsContainer(this->mpBoundaryConditions); // as the BCC can change
168  mpBidomainNeumannSurfaceTermAssembler->SetVectorToAssemble(this->mpLinearSystem->rGetRhsVector(), false/*don't zero vector!*/);
169  mpBidomainNeumannSurfaceTermAssembler->AssembleVector();
170 
171 
173  // apply correction term
175  if (mpBidomainCorrectionTermAssembler)
176  {
177  mpBidomainCorrectionTermAssembler->SetVectorToAssemble(this->mpLinearSystem->rGetRhsVector(), false/*don't zero vector!*/);
178  // don't need to set current solution
179  mpBidomainCorrectionTermAssembler->AssembleVector();
180  }
181 
182  this->mpLinearSystem->FinaliseRhsVector();
183 
184  this->mpBoundaryConditions->ApplyDirichletToLinearProblem(*(this->mpLinearSystem), computeMatrix);
185 
186  if (this->mBathSimulation)
187  {
188  this->mpLinearSystem->FinaliseLhsMatrix();
189  this->FinaliseForBath(computeMatrix,true);
190  }
191 
192  if (computeMatrix)
193  {
194  this->mpLinearSystem->FinaliseLhsMatrix();
195  }
196  this->mpLinearSystem->FinaliseRhsVector();
197 }
198 
199 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
201  bool bathSimulation,
203  BidomainTissue<SPACE_DIM>* pTissue,
205  : AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>(bathSimulation,pMesh,pTissue,pBoundaryConditions)
206 {
207  // Tell tissue there's no need to replicate ionic caches
208  pTissue->SetCacheReplication(false);
209  mVecForConstructingRhs = NULL;
210 
211  // create assembler
212  if (bathSimulation)
213  {
215  }
216  else
217  {
219  }
220 
221 
223 
225  {
228  //We are going to need those caches after all
229  pTissue->SetCacheReplication(true);
230  }
231  else
232  {
234  }
235 }
236 
237 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
239 {
240  delete mpBidomainAssembler;
241  delete mpBidomainNeumannSurfaceTermAssembler;
242 
243  if (mVecForConstructingRhs)
244  {
245  PetscTools::Destroy(mVecForConstructingRhs);
246  PetscTools::Destroy(mMassMatrix);
247  }
248 
249  if (mpBidomainCorrectionTermAssembler)
250  {
251  delete mpBidomainCorrectionTermAssembler;
252  }
253 }
254 
255 // Explicit instantiation
256 template class BidomainSolver<1,1>;
257 template class BidomainSolver<2,2>;
258 template class BidomainSolver<3,3>;
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
static bool IsRegionBath(HeartRegionType regionId)
BidomainAssembler< ELEMENT_DIM, SPACE_DIM > * mpBidomainAssembler
void InitialiseForSolve(Vec initialSolution)
DistributedVector CreateDistributedVector(Vec vec, bool readOnly=false)
BidomainNeumannSurfaceTermAssembler< ELEMENT_DIM, SPACE_DIM > * mpBidomainNeumannSurfaceTermAssembler
BidomainTissue< SPACE_DIM > * mpBidomainTissue
BidomainCorrectionTermAssembler< ELEMENT_DIM, SPACE_DIM > * mpBidomainCorrectionTermAssembler
bool GetUseStateVariableInterpolation() const
void SetupLinearSystem(Vec currentSolution, bool computeMatrix)
BidomainSolver(bool bathSimulation, AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, BidomainTissue< SPACE_DIM > *pTissue, BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, 2 > *pBoundaryConditions)
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
void SetCacheReplication(bool doCacheReplication)
static void Finalise(Mat matrix)
double GetCapacitance() const
double GetSurfaceAreaToVolumeRatio() const
void InitialiseForSolve(Vec initialSolution)
void SetMatrixToAssemble(Mat &rMatToAssemble, bool zeroMatrixBeforeAssembly=true)
static HeartConfig * Instance()