AbstractBidomainSolver.cpp

00001 
00002 /*
00003 
00004 Copyright (C) University of Oxford, 2005-2011
00005 
00006 University of Oxford means the Chancellor, Masters and Scholars of the
00007 University of Oxford, having an administrative office at Wellington
00008 Square, Oxford OX1 2JD, UK.
00009 
00010 This file is part of Chaste.
00011 
00012 Chaste is free software: you can redistribute it and/or modify it
00013 under the terms of the GNU Lesser General Public License as published
00014 by the Free Software Foundation, either version 2.1 of the License, or
00015 (at your option) any later version.
00016 
00017 Chaste is distributed in the hope that it will be useful, but WITHOUT
00018 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00019 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00020 License for more details. The offer of Chaste under the terms of the
00021 License is subject to the License being interpreted in accordance with
00022 English Law and subject to any action against the University of Oxford
00023 being under the jurisdiction of the English Courts.
00024 
00025 You should have received a copy of the GNU Lesser General Public License
00026 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00027 
00028 */
00029 
00030 
00031 #include "AbstractBidomainSolver.hpp"
00032 #include "TetrahedralMesh.hpp"
00033 #include "PetscMatTools.hpp"
00034 #include "PetscVecTools.hpp"
00035 
00036 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00037 void AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>::InitialiseForSolve(Vec initialSolution)
00038 {
00039     if (this->mpLinearSystem != NULL)
00040     {
00041         return;
00042     }
00043 
00044     // linear system created here
00045     AbstractDynamicLinearPdeSolver<ELEMENT_DIM,SPACE_DIM,2>::InitialiseForSolve(initialSolution);
00046 
00047     if (HeartConfig::Instance()->GetUseAbsoluteTolerance())
00048     {
00049 #ifdef TRACE_KSP
00050         if (PetscTools::AmMaster())
00051         {
00052             std::cout << "Using absolute tolerance: " << mpConfig->GetAbsoluteTolerance() <<"\n";
00053         }
00054 #endif
00055         this->mpLinearSystem->SetAbsoluteTolerance(mpConfig->GetAbsoluteTolerance());
00056     }
00057     else
00058     {
00059 #ifdef TRACE_KSP
00060         if (PetscTools::AmMaster())
00061         {
00062             std::cout << "Using relative tolerance: " << mpConfig->GetRelativeTolerance() <<"\n";
00063         }
00064 #endif
00065         this->mpLinearSystem->SetRelativeTolerance(mpConfig->GetRelativeTolerance());
00066     }
00067 
00068     this->mpLinearSystem->SetKspType(HeartConfig::Instance()->GetKSPSolver());
00069 
00071     if(std::string("twolevelsblockdiagonal") == std::string(HeartConfig::Instance()->GetKSPPreconditioner()))
00072     {
00074         TetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* p_mesh = dynamic_cast<TetrahedralMesh<ELEMENT_DIM,SPACE_DIM>*>(this->mpMesh);
00075         if (p_mesh && PetscTools::IsSequential())
00076         {
00077             boost::shared_ptr<std::vector<PetscInt> > p_bath_nodes(new std::vector<PetscInt>());
00078         
00079             for(unsigned node_index=0; node_index<this->mpMesh->GetNumNodes(); node_index++)
00080             {
00081                 if (HeartRegionCode::IsRegionBath( this->mpMesh->GetNode(node_index)->GetRegion() ))
00082                 {
00083                     p_bath_nodes->push_back(node_index);
00084                 }
00085             }
00086             
00087             this->mpLinearSystem->SetPcType(HeartConfig::Instance()->GetKSPPreconditioner(), p_bath_nodes);
00088         }
00089         else
00090         {
00091             TERMINATE("Two levels block diagonal only works with TetrahedralMesh and p=1");
00092         }
00093     }
00094     else
00095     { 
00096         this->mpLinearSystem->SetPcType(HeartConfig::Instance()->GetKSPPreconditioner());
00097     }
00098 
00099     if (mRowForAverageOfPhiZeroed==INT_MAX)
00100     {
00101         // not applying average(phi)=0 constraint, so matrix is symmetric
00102         this->mpLinearSystem->SetMatrixIsSymmetric(true);
00103     }
00104     else
00105     {
00106         // applying average(phi)=0 constraint, so matrix is not symmetric
00107         this->mpLinearSystem->SetMatrixIsSymmetric(false);
00108     }
00109 
00110     this->mpLinearSystem->SetUseFixedNumberIterations(HeartConfig::Instance()->GetUseFixedNumberIterationsLinearSolver(), HeartConfig::Instance()->GetEvaluateNumItsEveryNSolves());
00111 }
00112 
00113 
00114 
00115 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00116 void AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>::PrepareForSetupLinearSystem(Vec existingSolution)
00117 {
00118     double time = PdeSimulationTime::GetTime();
00119     double dt = PdeSimulationTime::GetPdeTimeStep();
00120     mpBidomainTissue->SolveCellSystems(existingSolution, time, time+dt);
00121 }
00122 
00123 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00124 Vec AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>::GenerateNullBasis() const
00125 {
00126     double sqrt_num_nodes = sqrt((double) this->mpMesh->GetNumNodes());
00127 
00128     Vec null_basis;
00129     DistributedVectorFactory* p_factory = this->mpMesh->GetDistributedVectorFactory();
00130     null_basis = p_factory->CreateVec(2);
00131 
00132     DistributedVector dist_null_basis = p_factory->CreateDistributedVector(null_basis);
00133     DistributedVector::Stripe null_basis_stripe_0(dist_null_basis,0);
00134     DistributedVector::Stripe null_basis_stripe_1(dist_null_basis,1);
00135     for (DistributedVector::Iterator index = dist_null_basis.Begin();
00136          index != dist_null_basis.End();
00137          ++index)
00138     {
00139         null_basis_stripe_0[index] = 0.0;
00140         null_basis_stripe_1[index] = 1.0/sqrt_num_nodes; // normalised vector
00141     }
00142     dist_null_basis.Restore();
00143 
00144     return null_basis;
00145 }
00146 
00147 
00148 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00149 void AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>::FinaliseLinearSystem(Vec existingSolution)
00150 {
00151     // If no dirichlet boundary conditions
00152     //  (i) Check compatibility condition to check we are solving
00153     //      a linear system that can be solved
00154     // Then either
00155     //  (a) If not setting average(phi)=0, we are solving a singular system,
00156     //      so set up a null space
00157     //  (b) Apply average(phi)=0 constraint by altering the last row, to
00158     //      get a non-singular system
00159     //
00160     if (!(GetBoundaryConditions()->HasDirichletBoundaryConditions()))
00161     {
00162         // first check compatibility condition
00163         CheckCompatibilityCondition();
00164 
00165         // Check whether applying average(phi_e)=0 constraint
00166         if (mRowForAverageOfPhiZeroed==INT_MAX)
00167         {
00168             // We're not using the constraint, hence use a null space
00169             if (!mNullSpaceCreated)
00170             {
00171                 // No null space set up, so create one and pass it to the linear system
00172                 Vec null_basis[] = {GenerateNullBasis()};
00173 
00174                 this->mpLinearSystem->SetNullBasis(null_basis, 1);
00175 
00176                 VecDestroy(null_basis[0]);
00177                 mNullSpaceCreated = true;
00178             }
00179         }
00180         else
00181         {
00182             // Using the average(phi_e)=0 constraint
00183 
00184             // CG (default solver) won't work since the system isn't symmetric anymore. Switch to GMRES
00185             this->mpLinearSystem->SetKspType("gmres"); // Switches the solver
00186             mpConfig->SetKSPSolver("gmres"); // Makes sure this change will be reflected in the XML file written to disk at the end of the simulation.
00187 
00188             // Set average phi_e to zero
00189             unsigned matrix_size = this->mpLinearSystem->GetSize();
00190             if (!this->mMatrixIsAssembled)
00191             {
00192 
00193                 // Set the mRowForAverageOfPhiZeroed-th matrix row to 0 1 0 1 ...
00194                 std::vector<unsigned> row_for_average;
00195                 row_for_average.push_back(mRowForAverageOfPhiZeroed);
00196                 this->mpLinearSystem->ZeroMatrixRowsWithValueOnDiagonal(row_for_average, 0.0);
00197                 for (unsigned col_index=0; col_index<matrix_size; col_index++)
00198                 {
00199                     if (col_index%2 == 1)
00200                     {
00201                         this->mpLinearSystem->SetMatrixElement(mRowForAverageOfPhiZeroed, col_index, 1);
00202                     }
00203 
00204                 }
00205                 this->mpLinearSystem->AssembleFinalLhsMatrix();
00206 
00207             }
00208             // Set the mRowForAverageOfPhiZeroed-th rhs vector row to 0
00209             this->mpLinearSystem->SetRhsVectorElement(mRowForAverageOfPhiZeroed, 0);
00210 
00211             this->mpLinearSystem->AssembleRhsVector();
00212         }
00213     }
00214 }
00215 
00216 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00217 void AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>::CheckCompatibilityCondition()
00218 {
00219     if (GetBoundaryConditions()->HasDirichletBoundaryConditions() || this->mRowForAverageOfPhiZeroed!=INT_MAX)
00220     {
00221         // not a singular system, no compatibility condition
00222         return;
00223     }
00224 
00225 #ifndef NDEBUG   
00226     DistributedVector distributed_rhs=this->mpMesh->GetDistributedVectorFactory()->CreateDistributedVector(this->mpLinearSystem->rGetRhsVector());
00227     DistributedVector::Stripe distributed_rhs_phi_entries(distributed_rhs,1);
00228 
00229     double local_sum=0;
00230     for (DistributedVector::Iterator index = distributed_rhs.Begin();
00231          index!= distributed_rhs.End();
00232          ++index)
00233     {
00234         local_sum += distributed_rhs_phi_entries[index];
00235     }
00236     
00237     double global_sum;
00238     int mpi_ret = MPI_Allreduce(&local_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, PETSC_COMM_WORLD);
00239     assert(mpi_ret == MPI_SUCCESS);
00240      
00241     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)
00242     {
00243         #define COVERAGE_IGNORE
00244         // shouldn't ever reach this line but useful to have the error printed out if you do
00245         //std::cout << "Sum of b_{2i+1} = " << global_sum << " (should be zero for compatibility)\n";
00246         EXCEPTION("Linear system does not satisfy compatibility constraint!");
00247         #undef COVERAGE_IGNORE
00248     }
00249 #endif
00250 }
00251 
00252 
00253 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00254 AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>::AbstractBidomainSolver(
00255             bool bathSimulation,
00256             AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00257             BidomainTissue<SPACE_DIM>* pTissue,
00258             BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,2>* pBoundaryConditions,
00259             unsigned numQuadPoints)
00260     : AbstractDynamicLinearPdeSolver<ELEMENT_DIM,SPACE_DIM,2>(pMesh),
00261       mBathSimulation(bathSimulation),
00262       mpBidomainTissue(pTissue),
00263       mpBoundaryConditions(pBoundaryConditions),
00264       mNumQuadPoints(numQuadPoints)
00265 {
00266     assert(pTissue != NULL);
00267     assert(pBoundaryConditions != NULL);
00268 
00269     mNullSpaceCreated = false;
00270 
00271     // important!
00272     this->mMatrixIsConstant = true;
00273 
00274     mRowForAverageOfPhiZeroed = INT_MAX; //this->mpLinearSystem->GetSize() - 1;
00275     mpConfig = HeartConfig::Instance();
00276 
00277     mpBidomainAssembler = NULL; // can't initialise until know what dt is
00278 }
00279 
00280 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00281 AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>::~AbstractBidomainSolver()
00282 {
00283     if(mpBidomainAssembler)
00284     {
00285         delete mpBidomainAssembler;
00286     }
00287 }
00288 
00289 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00290 void AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>::SetFixedExtracellularPotentialNodes(
00291             std::vector<unsigned> fixedExtracellularPotentialNodes)
00292 {
00293     for (unsigned i=0; i<fixedExtracellularPotentialNodes.size(); i++)
00294     {
00295         if (fixedExtracellularPotentialNodes[i] >= this->mpMesh->GetNumNodes() )
00296         {
00297             EXCEPTION("Fixed node number must be less than total number nodes");
00298         }
00299     }
00300 
00301     mFixedExtracellularPotentialNodes = fixedExtracellularPotentialNodes;
00302 
00303     // We will need to recalculate this when HasDirichletBoundaryConditions() is called.
00304     GetBoundaryConditions()->ResetDirichletCommunication();
00305 
00306     for (unsigned i=0; i<mFixedExtracellularPotentialNodes.size(); i++)
00307     {
00308         if (this->mpMesh->GetDistributedVectorFactory()->IsGlobalIndexLocal(mFixedExtracellularPotentialNodes[i]))
00309         {
00310             ConstBoundaryCondition<SPACE_DIM>* p_boundary_condition
00311                  = new ConstBoundaryCondition<SPACE_DIM>(0.0);
00312 
00313             //Throws if node is not owned locally
00314             Node<SPACE_DIM>* p_node = this->mpMesh->GetNode(mFixedExtracellularPotentialNodes[i]);
00315 
00316             GetBoundaryConditions()->AddDirichletBoundaryCondition(p_node, p_boundary_condition, 1);
00317 
00318         }
00319     }
00320 }
00321 
00322 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00323 void AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>::SetRowForAverageOfPhiZeroed(unsigned row)
00324 {
00325     // Row should be odd in C++-like indexing
00326     if (row%2 == 0)
00327     {
00328         EXCEPTION("Row for applying the constraint 'Average of phi_e = zero' should be odd in C++ like indexing");
00329     }
00330 
00331     mRowForAverageOfPhiZeroed = row;
00332 }
00333 
00334 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00335 void AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>::FinaliseForBath(bool computeMatrix, bool computeVector)
00336 {
00337     assert(mBathSimulation);
00338 
00339     PetscTruth is_matrix_assembled;
00340     MatAssembled(this->mpLinearSystem->GetLhsMatrix(), &is_matrix_assembled);
00341     assert(is_matrix_assembled == PETSC_TRUE);
00342 
00343     /*
00344      *  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.
00345      *  When assembling a bath element you get a matrix subblock that looks like (2D example):
00346      *
00347      *     Vm   0 0 0 0 0 0
00348      *     Vm   0 0 0 0 0 0
00349      *     Vm   0 0 0 0 0 0
00350      *     Phie 0 0 0 x x x
00351      *     Phie 0 0 0 x x x  -> the x subblock is assembled from div(grad_phi) = 0
00352      *     Phie 0 0 0 x x x
00353      *
00354      *  Therefore, all the Vm entries of this node are already 0.
00355      *
00356      *  Explicitly checking it in non-production builds.
00357      */
00358 #ifndef NDEBUG
00359     if(computeMatrix)
00360     {
00361         for (typename AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>::NodeIterator iter=this->mpMesh->GetNodeIteratorBegin();
00362              iter != this->mpMesh->GetNodeIteratorEnd();
00363              ++iter)
00364         {
00365             if (HeartRegionCode::IsRegionBath( (*iter).GetRegion() ))
00366             {
00367                 int num_equation = 2*iter->GetIndex(); // assumes Vm and Phie are interleaved
00368 
00369                 PetscInt local_lo, local_hi;
00370                 this->mpLinearSystem->GetOwnershipRange(local_lo, local_hi);
00371 
00372                 // If this processor owns i-th row, check it.
00373                 if ((local_lo <= (int)num_equation) && ((int)num_equation < local_hi))
00374                 {
00375                     for (unsigned column=0; column < this->mpLinearSystem->GetSize(); column++)
00376                     {
00377                         assert(this->mpLinearSystem->GetMatrixElement(num_equation, column)==0.0);
00378                     }
00379                 }
00380 
00381                 // Check the local entries of the i-th column
00382                 for (int row=local_lo; row<local_hi; row++)
00383                 {
00384                     assert(this->mpLinearSystem->GetMatrixElement(row, num_equation)==0);
00385                 }
00386             }
00387         }
00388     }
00389 #endif
00390 
00391     /*
00392      *  These two loops are decoupled because interleaving calls to GetMatrixElement and MatSetValue
00393      *  require reassembling the matrix before GetMatrixElement which generates massive communication
00394      *  overhead for large models and/or large core counts.
00395      */
00396 
00397     for (typename AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>::NodeIterator iter=this->mpMesh->GetNodeIteratorBegin();
00398          iter != this->mpMesh->GetNodeIteratorEnd();
00399          ++iter)
00400     {
00401         if (HeartRegionCode::IsRegionBath((*iter).GetRegion() ))
00402         {
00403             PetscInt index = 2*iter->GetIndex(); // assumes Vm and Phie are interleaved
00404 
00405             if(computeMatrix)
00406             {
00407                 // put 1.0 on the diagonal
00408                 PetscMatTools::SetElement(this->mpLinearSystem->rGetLhsMatrix(), index, index, 1.0);
00409             }
00410 
00411             if(computeVector)
00412             {
00413                 // zero rhs vector entry
00414                 PetscVecTools::SetElement(this->mpLinearSystem->rGetRhsVector(), index, 0.0);
00415             }
00416         }
00417     }
00418 }
00419 
00421 // explicit instantiation
00423 
00424 template class AbstractBidomainSolver<1,1>;
00425 template class AbstractBidomainSolver<2,2>;
00426 template class AbstractBidomainSolver<3,3>;

Generated on Tue May 31 14:31:41 2011 for Chaste by  doxygen 1.5.5