Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
AbstractBidomainSolver.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
37#include "AbstractBidomainSolver.hpp"
38#include "TetrahedralMesh.hpp"
39#include "PetscMatTools.hpp"
40#include "PetscVecTools.hpp"
41
42template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
44{
45 // The base class method that calls this function will only call it with a null linear system
46 assert(this->mpLinearSystem == NULL);
47
48 // linear system created here
50
51 if (HeartConfig::Instance()->GetUseAbsoluteTolerance())
52 {
53#ifdef TRACE_KSP
55 {
56 std::cout << "Using absolute tolerance: " << mpConfig->GetAbsoluteTolerance() << "\n";
57 }
58#endif
59 this->mpLinearSystem->SetAbsoluteTolerance(mpConfig->GetAbsoluteTolerance());
60 }
61 else
62 {
63#ifdef TRACE_KSP
65 {
66 std::cout << "Using relative tolerance: " << mpConfig->GetRelativeTolerance() << "\n";
67 }
68#endif
69 this->mpLinearSystem->SetRelativeTolerance(mpConfig->GetRelativeTolerance());
70 }
71
72 this->mpLinearSystem->SetKspType(HeartConfig::Instance()->GetKSPSolver());
73
74 // Note that twolevelblockdiagonal was never finished. It was a preconditioner specific to the Parabolic-Parabolic formulation of Bidomain
75 // Two levels block diagonal only worked in serial with TetrahedralMesh.
76 assert(std::string(HeartConfig::Instance()->GetKSPPreconditioner()) != std::string("twolevelsblockdiagonal"));
77 this->mpLinearSystem->SetPcType(HeartConfig::Instance()->GetKSPPreconditioner());
78
79 if (mRowForAverageOfPhiZeroed == INT_MAX)
80 {
81 // not applying average(phi)=0 constraint, so matrix is symmetric
82 this->mpLinearSystem->SetMatrixIsSymmetric(true);
83 }
84 else
85 {
86 //Turn off preallocation so that we can have one dense row in the matrix.
87 PetscMatTools::TurnOffVariableAllocationError(this->mpLinearSystem->rGetLhsMatrix());
88
89 // applying average(phi)=0 constraint, so matrix is not symmetric
90 this->mpLinearSystem->SetMatrixIsSymmetric(false);
91 }
92
93 this->mpLinearSystem->SetUseFixedNumberIterations(
94 HeartConfig::Instance()->GetUseFixedNumberIterationsLinearSolver(),
95 HeartConfig::Instance()->GetEvaluateNumItsEveryNSolves());
96}
97
98template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
100{
101 mpBidomainTissue->SolveCellSystems(existingSolution, PdeSimulationTime::GetTime(), PdeSimulationTime::GetNextTime());
102}
103
104template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
106{
107 double sqrt_num_nodes = sqrt((double) this->mpMesh->GetNumNodes());
108
109 Vec null_basis;
110 DistributedVectorFactory* p_factory = this->mpMesh->GetDistributedVectorFactory();
111 null_basis = p_factory->CreateVec(2);
112
113 DistributedVector dist_null_basis = p_factory->CreateDistributedVector(null_basis);
114 DistributedVector::Stripe null_basis_stripe_0(dist_null_basis,0);
115 DistributedVector::Stripe null_basis_stripe_1(dist_null_basis,1);
116 for (DistributedVector::Iterator index = dist_null_basis.Begin();
117 index != dist_null_basis.End();
118 ++index)
119 {
120 null_basis_stripe_0[index] = 0.0;
121 null_basis_stripe_1[index] = 1.0/sqrt_num_nodes; // normalised vector
122 }
123 dist_null_basis.Restore();
124
125 return null_basis;
126}
127
128template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
130{
131 // If no dirichlet boundary conditions
132 // (i) Check compatibility condition to check we are solving
133 // a linear system that can be solved
134 // Then either
135 // (a) If not setting average(phi)=0, we are solving a singular system,
136 // so set up a null space
137 // (b) Apply average(phi)=0 constraint by altering the last row, to
138 // get a non-singular system
139 //
140 if (!(GetBoundaryConditions()->HasDirichletBoundaryConditions()))
141 {
142 // first check compatibility condition
143 CheckCompatibilityCondition();
144
145 // Check whether applying average(phi_e)=0 constraint
146 if (mRowForAverageOfPhiZeroed==INT_MAX)
147 {
148 // We're not using the constraint, hence use a null space
149 if (!mNullSpaceCreated)
150 {
151 // No null space set up, so create one and pass it to the linear system
152 Vec null_basis[] = {GenerateNullBasis()};
153
154 this->mpLinearSystem->SetNullBasis(null_basis, 1);
155
156 PetscTools::Destroy(null_basis[0]);
157 mNullSpaceCreated = true;
158 }
160 else
161 {
162 // Using the average(phi_e)=0 constraint
163
164 // CG (default solver) won't work since the system isn't symmetric anymore. Switch to GMRES
165 this->mpLinearSystem->SetKspType("gmres"); // Switches the solver
166 mpConfig->SetKSPSolver("gmres", true); // Makes sure this change will be reflected in the XML file written to disk at the end of the simulation.
167 //(If the user doesn't have gmres then the "true" warns the user about the switch)
168
169 // Set average phi_e to zero
170 unsigned matrix_size = this->mpLinearSystem->GetSize();
171 if (!this->mMatrixIsAssembled)
172 {
173
174 // Set the mRowForAverageOfPhiZeroed-th matrix row to 0 1 0 1 ...
175 std::vector<unsigned> row_for_average;
176 row_for_average.push_back(mRowForAverageOfPhiZeroed);
177 this->mpLinearSystem->ZeroMatrixRowsWithValueOnDiagonal(row_for_average, 0.0);
178 for (unsigned col_index=0; col_index<matrix_size; col_index++)
180 if (col_index%2 == 1)
181 {
182 this->mpLinearSystem->SetMatrixElement(mRowForAverageOfPhiZeroed, col_index, 1);
183 }
184
185 }
186 this->mpLinearSystem->FinaliseLhsMatrix();
188 }
189 // Set the mRowForAverageOfPhiZeroed-th rhs vector row to 0
190 this->mpLinearSystem->SetRhsVectorElement(mRowForAverageOfPhiZeroed, 0);
191
192 this->mpLinearSystem->FinaliseRhsVector();
193 }
194 }
195}
196
197template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
199{
200 if (GetBoundaryConditions()->HasDirichletBoundaryConditions() || this->mRowForAverageOfPhiZeroed!=INT_MAX)
201 {
202 // not a singular system, no compatibility condition
203 return;
204 }
205
206#ifndef NDEBUG
207 DistributedVector distributed_rhs=this->mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(this->mpLinearSystem->rGetRhsVector());
208 DistributedVector::Stripe distributed_rhs_phi_entries(distributed_rhs,1);
209
210 double local_sum=0;
211 for (DistributedVector::Iterator index = distributed_rhs.Begin();
212 index!= distributed_rhs.End();
213 ++index)
214 {
215 local_sum += distributed_rhs_phi_entries[index];
216 }
217
218 double global_sum;
219 int mpi_ret = MPI_Allreduce(&local_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
220 assert(mpi_ret == MPI_SUCCESS);
221
222 if (fabs(global_sum)>1e-6) // magic number! sum should really be a sum of zeros and exactly zero though anyway (or a-a+b-b+c-c.. etc in the case of electrodes)
223 {
224 // LCOV_EXCL_START
225 // shouldn't ever reach this line but useful to have the error printed out if you do
226 //std::cout << "Sum of b_{2i+1} = " << global_sum << " (should be zero for compatibility)\n";
227 EXCEPTION("Linear system does not satisfy compatibility constraint!");
228 // LCOV_EXCL_STOP
229 }
230#endif
231}
232
233template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
235 bool bathSimulation,
239 : AbstractDynamicLinearPdeSolver<ELEMENT_DIM,SPACE_DIM,2>(pMesh),
240 mBathSimulation(bathSimulation),
241 mpBidomainTissue(pTissue),
242 mpBoundaryConditions(pBoundaryConditions)
243{
244 assert(pTissue != NULL);
245 assert(pBoundaryConditions != NULL);
246
247 mNullSpaceCreated = false;
248
249 // important!
250 this->mMatrixIsConstant = true;
251
252 mRowForAverageOfPhiZeroed = INT_MAX; //this->mpLinearSystem->GetSize() - 1;
254}
255
256template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
260
261template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
263 std::vector<unsigned> fixedExtracellularPotentialNodes)
264{
265 for (unsigned i=0; i<fixedExtracellularPotentialNodes.size(); i++)
266 {
267 if (fixedExtracellularPotentialNodes[i] >= this->mpMesh->GetNumNodes() )
268 {
269 EXCEPTION("Fixed node number must be less than total number nodes");
270 }
271 }
272
273 mFixedExtracellularPotentialNodes = fixedExtracellularPotentialNodes;
274
275 // We will need to recalculate this when HasDirichletBoundaryConditions() is called.
276 GetBoundaryConditions()->ResetDirichletCommunication();
277
278 for (unsigned i=0; i<mFixedExtracellularPotentialNodes.size(); i++)
279 {
280 if (this->mpMesh->GetDistributedVectorFactory()->IsGlobalIndexLocal(mFixedExtracellularPotentialNodes[i]))
281 {
282 ConstBoundaryCondition<SPACE_DIM>* p_boundary_condition
284
285 //Throws if node is not owned locally
286 Node<SPACE_DIM>* p_node = this->mpMesh->GetNode(mFixedExtracellularPotentialNodes[i]);
287
288 GetBoundaryConditions()->AddDirichletBoundaryCondition(p_node, p_boundary_condition, 1);
289
290 }
291 }
292}
293
294template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
296{
297 // Row should be odd in C++-like indexing
298 if (row%2 == 0)
299 {
300 EXCEPTION("Row for applying the constraint 'Average of phi_e = zero' should be odd in C++ like indexing");
301 }
302
303 mRowForAverageOfPhiZeroed = row;
304}
305
306template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
307void AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>::FinaliseForBath(bool computeMatrix, bool computeVector)
308{
309 assert(mBathSimulation);
310 PetscBool is_matrix_assembled;
311 MatAssembled(this->mpLinearSystem->GetLhsMatrix(), &is_matrix_assembled);
312 assert(is_matrix_assembled == PETSC_TRUE);
313
314 /*
315 * Before revision 6516, we used to zero out i-th row and column here. It seems to be redundant because they are already zero after assembly.
316 * When assembling a bath element you get a matrix subblock that looks like (2D example):
317 *
318 * Vm 0 0 0 0 0 0
319 * Vm 0 0 0 0 0 0
320 * Vm 0 0 0 0 0 0
321 * Phie 0 0 0 x x x
322 * Phie 0 0 0 x x x -> the x subblock is assembled from div(grad_phi) = 0
323 * Phie 0 0 0 x x x
324 *
325 * Therefore, all the Vm entries of this node are already 0.
326 *
327 * Explicitly checking it in non-production builds.
328 */
329#ifndef NDEBUG
330 if (computeMatrix)
331 {
332 for (typename AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>::NodeIterator iter=this->mpMesh->GetNodeIteratorBegin();
333 iter != this->mpMesh->GetNodeIteratorEnd();
334 ++iter)
335 {
336 if (HeartRegionCode::IsRegionBath( (*iter).GetRegion() ))
337 {
338 int num_equation = 2*iter->GetIndex(); // assumes Vm and Phie are interleaved
339
340 PetscInt local_lo, local_hi;
341 this->mpLinearSystem->GetOwnershipRange(local_lo, local_hi);
342
343 // If this processor owns i-th row, check it.
344 if ((local_lo <= (int)num_equation) && ((int)num_equation < local_hi))
345 {
346 for (unsigned column=0; column < this->mpLinearSystem->GetSize(); column++)
347 {
348 assert(this->mpLinearSystem->GetMatrixElement(num_equation, column)==0.0);
349 }
350 }
351
352 // Check the local entries of the i-th column
353 for (int row=local_lo; row<local_hi; row++)
354 {
355 assert(this->mpLinearSystem->GetMatrixElement(row, num_equation)==0);
356 }
357 }
358 }
359 }
360#endif
361
362 /*
363 * These two loops are decoupled because interleaving calls to GetMatrixElement and MatSetValue
364 * require reassembling the matrix before GetMatrixElement which generates massive communication
365 * overhead for large models and/or large core counts.
366 */
367
368 for (typename AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>::NodeIterator iter=this->mpMesh->GetNodeIteratorBegin();
369 iter != this->mpMesh->GetNodeIteratorEnd();
370 ++iter)
371 {
372 if (HeartRegionCode::IsRegionBath((*iter).GetRegion() ))
373 {
374 PetscInt index = 2*iter->GetIndex(); // assumes Vm and Phie are interleaved
375
376 if (computeMatrix)
377 {
378 // put 1.0 on the diagonal
379 PetscMatTools::SetElement(this->mpLinearSystem->rGetLhsMatrix(), index, index, 1.0);
380 }
381
382 if (computeVector)
383 {
384 // zero rhs vector entry
385 PetscVecTools::SetElement(this->mpLinearSystem->rGetRhsVector(), index, 0.0);
386 }
387 }
388 }
389}
390
391// Explicit instantiation
392template class AbstractBidomainSolver<1,1>;
393template class AbstractBidomainSolver<2,2>;
394template class AbstractBidomainSolver<3,3>;
#define EXCEPTION(message)
void InitialiseForSolve(Vec initialSolution)
virtual Vec GenerateNullBasis() const
void SetRowForAverageOfPhiZeroed(unsigned rowMeanPhiEZero)
void FinaliseForBath(bool computeMatrix, bool computeVector)
AbstractBidomainSolver(bool bathSimulation, AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, BidomainTissue< SPACE_DIM > *pTissue, BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, 2 > *pBoundaryConditions)
void SetFixedExtracellularPotentialNodes(std::vector< unsigned > fixedExtracellularPotentialNodes)
virtual void FinaliseLinearSystem(Vec existingSolution)
void PrepareForSetupLinearSystem(Vec existingSolution)
virtual void InitialiseForSolve(Vec initialSolution=nullptr)
DistributedVector CreateDistributedVector(Vec vec, bool readOnly=false)
static HeartConfig * Instance()
static bool IsRegionBath(HeartRegionType regionId)
Definition Node.hpp:59
static double GetTime()
static double GetNextTime()
static void SetElement(Mat matrix, PetscInt row, PetscInt col, double value)
static void TurnOffVariableAllocationError(Mat matrix)
static void Destroy(Vec &rVec)
static bool AmMaster()
static void SetElement(Vec vector, PetscInt row, double value)