Chaste Commit::baa90ac2819b962188b7562f2326be23c47859a7
AbstractExtendedBidomainSolver.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 "AbstractExtendedBidomainSolver.hpp"
37
38template <unsigned ELEMENT_DIM, unsigned SPACE_DIM>
40{
41 // The base class method that calls this function will only call it with a null linear system
42 assert(this->mpLinearSystem == NULL);
43
44 // linear system created here
46
47 if (HeartConfig::Instance()->GetUseAbsoluteTolerance())
48 {
49#ifdef TRACE_KSP
50 std::cout << "Using absolute tolerance: " << mpConfig->GetAbsoluteTolerance() << "\n";
51#endif
52 this->mpLinearSystem->SetAbsoluteTolerance(mpConfig->GetAbsoluteTolerance());
53 }
54 else
55 {
56#ifdef TRACE_KSP
57 std::cout << "Using relative tolerance: " << mpConfig->GetRelativeTolerance() << "\n";
58#endif
59 this->mpLinearSystem->SetRelativeTolerance(mpConfig->GetRelativeTolerance());
60 }
61
62 this->mpLinearSystem->SetKspType(HeartConfig::Instance()->GetKSPSolver());
63 this->mpLinearSystem->SetPcType(HeartConfig::Instance()->GetKSPPreconditioner());
64
65 if (mRowForAverageOfPhiZeroed == INT_MAX)
66 {
67 // not applying average(phi)=0 constraint, so matrix is symmetric
68 this->mpLinearSystem->SetMatrixIsSymmetric(true);
69 }
70 else
71 {
72 //Turn off preallocation so that we can have one dense row in the matrix.
73 PetscMatTools::TurnOffVariableAllocationError(this->mpLinearSystem->rGetLhsMatrix());
74 // applying average(phi)=0 constraint, so matrix is not symmetric
75 this->mpLinearSystem->SetMatrixIsSymmetric(false);
76 }
77}
78
79template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
81{
82 mpExtendedBidomainTissue->SolveCellSystems(existingSolution, PdeSimulationTime::GetTime(), PdeSimulationTime::GetNextTime());
83}
84
85template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
87{
88 double sqrrt_num_nodes = sqrt( this->mpMesh->GetNumNodes());
89 double normalised_basis_value = 1.0 / sqrrt_num_nodes;
90
91 Vec nullbasis;
92 DistributedVectorFactory* p_factory = this->mpMesh->GetDistributedVectorFactory();
93 nullbasis=p_factory->CreateVec(3);
94 DistributedVector dist_null_basis = p_factory->CreateDistributedVector(nullbasis);
95 DistributedVector::Stripe null_basis_stripe_0(dist_null_basis,0);
96 DistributedVector::Stripe null_basis_stripe_1(dist_null_basis,1);
97 DistributedVector::Stripe null_basis_stripe_2(dist_null_basis,2);
98 for (DistributedVector::Iterator index = dist_null_basis.Begin();
99 index != dist_null_basis.End();
100 ++index)
101 {
102 null_basis_stripe_0[index] = 0.0;
103 null_basis_stripe_1[index] = 0.0;
104 null_basis_stripe_2[index] = normalised_basis_value;
105 }
106 dist_null_basis.Restore();
107 return nullbasis;
108}
109
110template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
112{
113 if (!(this->mpBoundaryConditions->HasDirichletBoundaryConditions()))
114 {
115 // We're not pinning any nodes.
116 if (mRowForAverageOfPhiZeroed==INT_MAX)
117 {
118 // We're not using the 'Average phi_e = 0' method, hence use a null space
119 if (!mNullSpaceCreated)
120 {
121 // No null space set up, so create one and pass it to the linear system
122 Vec nullbasis[] = {GenerateNullBasis()};
124 this->mpLinearSystem->SetNullBasis(nullbasis, 1);
125
126 PetscTools::Destroy(nullbasis[0]);
127 mNullSpaceCreated = true;
128 }
129 }
130 else // mRowForAverageOfPhiZeroed!=INT_MAX, i.e. we're using the 'Average phi_e = 0' method
131 {
132 // CG (default solver) won't work since the system isn't symmetric anymore. Switch to GMRES
133 this->mpLinearSystem->SetKspType("gmres"); // Switches the solver
134 mpConfig->SetKSPSolver("gmres", true); // Makes sure this change will be reflected in the XML file written to disk at the end of the simulation.
135 //(If the user doesn't have gmres then the "true" warns the user about the switch)
136
137 // Set average phi_e to zero
138 unsigned matrix_size = this->mpLinearSystem->GetSize();
139 if (!this->mMatrixIsAssembled)
140 {
141
142 // Set the mRowForAverageOfPhiZeroed-th matrix row to 0 0 1 0 0 1 ...
143 std::vector<unsigned> row_for_average;
144 row_for_average.push_back(mRowForAverageOfPhiZeroed);
145 this->mpLinearSystem->ZeroMatrixRowsWithValueOnDiagonal(row_for_average, 0.0);
146 for (unsigned col_index=0; col_index<matrix_size; col_index++)
147 {
148 if (((col_index-2)%3 == 0) && (col_index>1))//phi_e column indices are 2,5,8,11 etc....
149 {
150 this->mpLinearSystem->SetMatrixElement(mRowForAverageOfPhiZeroed, col_index, 1);
151 }
152
153 }
154 this->mpLinearSystem->FinaliseLhsMatrix();
155 }
156 // Set the mRowForAverageOfPhiZeroed-th rhs vector row to 0
157 this->mpLinearSystem->SetRhsVectorElement(mRowForAverageOfPhiZeroed, 0);
158
159 this->mpLinearSystem->FinaliseRhsVector();
160 }
161 }
162 CheckCompatibilityCondition();
163}
164
165template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
167{
168 if (this->mpBoundaryConditions->HasDirichletBoundaryConditions() || mRowForAverageOfPhiZeroed!=INT_MAX )
169 {
170 // not a singular system, no compability condition
171 return;
172 }
173
174#ifndef NDEBUG
175 DistributedVector distributed_rhs=this->mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(this->mpLinearSystem->rGetRhsVector());
176 DistributedVector::Stripe distributed_rhs_phi_entries(distributed_rhs,2); // stripe number 2 -> phi_e
177
178 double local_sum=0;
179 for (DistributedVector::Iterator index = distributed_rhs.Begin();
180 index!= distributed_rhs.End();
181 ++index)
182 {
183 local_sum += distributed_rhs_phi_entries[index];
184 }
185
186 double global_sum;
187 int mpi_ret = MPI_Allreduce(&local_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
188 assert(mpi_ret == MPI_SUCCESS);
189
190 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)
191 {
192 // LCOV_EXCL_START
193 // shouldn't ever reach this line but useful to have the error printed out if you do
194 std::cout << "Sum of b_{every 3 items} = " << global_sum << " (should be zero for compatibility)\n";
195 EXCEPTION("Linear system does not satisfy compatibility constraint!");
196 // LCOV_EXCL_STOP
197 }
198#endif
199}
200
201template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
203 bool bathSimulation,
207 : AbstractDynamicLinearPdeSolver<ELEMENT_DIM,SPACE_DIM,3>(pMesh),
208 mBathSimulation(bathSimulation),
209 mpExtendedBidomainTissue(pTissue),
210 mpBoundaryConditions(pBcc)
211{
212 assert(pTissue != NULL);
213 assert(pBcc != NULL);
214
215 mNullSpaceCreated = false;
216
217 // important!
218 this->mMatrixIsConstant = true;
219
220 mRowForAverageOfPhiZeroed = INT_MAX; //this->mpLinearSystem->GetSize() - 1;
222
223 mpExtendedBidomainAssembler = NULL; // can't initialise until know what dt is
224}
225
226template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
230
231template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
233 std::vector<unsigned> fixedExtracellularPotentialNodes)
234{
235 for (unsigned i=0; i<fixedExtracellularPotentialNodes.size(); i++)
236 {
237 if (fixedExtracellularPotentialNodes[i] >= this->mpMesh->GetNumNodes() )
238 {
239 EXCEPTION("Fixed node number must be less than total number nodes");
240 }
241 }
242
243 mFixedExtracellularPotentialNodes = fixedExtracellularPotentialNodes;
244
245 // We will need to recalculate this when HasDirichletBoundaryConditions() is called.
246 this->mpBoundaryConditions->ResetDirichletCommunication();
247
248 for (unsigned i=0; i<mFixedExtracellularPotentialNodes.size(); i++)
249 {
250 if (this->mpMesh->GetDistributedVectorFactory()->IsGlobalIndexLocal(mFixedExtracellularPotentialNodes[i]))
251 {
253
254 //Throws if node is not owned locally
255 Node<SPACE_DIM>* p_node = this->mpMesh->GetNode(mFixedExtracellularPotentialNodes[i]);
256
257 //the "false" passed in tells not to check that it is a boundary node (since we might have grounded electrodes within the tissue)
258 GetBoundaryConditions()->AddDirichletBoundaryCondition(p_node, p_boundary_condition, 2, false);
259 }
260 }
261}
262
263template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
265{
266 // Row should be every 3 entries, starting from zero...
267 if (((row-2)%3 != 0) && row!=INT_MAX)
268 {
269 EXCEPTION("Row for applying the constraint 'Average of phi_e = zero' should be every 3 rows");
270 }
271
272 mRowForAverageOfPhiZeroed = row;
273}
274
275// Explicit instantiation
#define EXCEPTION(message)
void SetFixedExtracellularPotentialNodes(std::vector< unsigned > fixedExtracellularPotentialNodes)
ExtendedBidomainAssembler< ELEMENT_DIM, SPACE_DIM > * mpExtendedBidomainAssembler
virtual void FinaliseLinearSystem(Vec existingSolution)
void PrepareForSetupLinearSystem(Vec existingSolution)
void SetRowForAverageOfPhiZeroed(unsigned rowMeanPhiEZero)
AbstractExtendedBidomainSolver(bool bathSimulation, AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, ExtendedBidomainTissue< SPACE_DIM > *pTissue, BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, 3 > *pBcc)
virtual void InitialiseForSolve(Vec initialSolution=nullptr)
DistributedVector CreateDistributedVector(Vec vec, bool readOnly=false)
static HeartConfig * Instance()
Definition Node.hpp:59
static double GetTime()
static double GetNextTime()
static void TurnOffVariableAllocationError(Mat matrix)
static void Destroy(Vec &rVec)