MatrixBasedBidomainSolver.cpp

00001 /*
00002 
00003 Copyright (C) University of Oxford, 2005-2011
00004 
00005 University of Oxford means the Chancellor, Masters and Scholars of the
00006 University of Oxford, having an administrative office at Wellington
00007 Square, Oxford OX1 2JD, UK.
00008 
00009 This file is part of Chaste.
00010 
00011 Chaste is free software: you can redistribute it and/or modify it
00012 under the terms of the GNU Lesser General Public License as published
00013 by the Free Software Foundation, either version 2.1 of the License, or
00014 (at your option) any later version.
00015 
00016 Chaste is distributed in the hope that it will be useful, but WITHOUT
00017 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
00018 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
00019 License for more details. The offer of Chaste under the terms of the
00020 License is subject to the License being interpreted in accordance with
00021 English Law and subject to any action against the University of Oxford
00022 being under the jurisdiction of the English Courts.
00023 
00024 You should have received a copy of the GNU Lesser General Public License
00025 along with Chaste. If not, see <http://www.gnu.org/licenses/>.
00026 
00027 */
00028 
00029 
00030 #include "MatrixBasedBidomainSolver.hpp"
00031 #include "BidomainAssembler.hpp"
00032 #include "BidomainWithBathAssembler.hpp"
00033 #include "PetscMatTools.hpp"
00034 
00035 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00036 void MatrixBasedBidomainSolver<ELEMENT_DIM,SPACE_DIM>::InitialiseForSolve(Vec initialSolution)
00037 {
00038     if (this->mpLinearSystem != NULL)
00039     {
00040         return;
00041     }
00042     AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>::InitialiseForSolve(initialSolution);
00043 
00044     // initialise matrix-based RHS vector and matrix, and use the linear
00045     // system rhs as a template
00046     Vec& r_template = this->mpLinearSystem->rGetRhsVector();
00047     VecDuplicate(r_template, &mVecForConstructingRhs);
00048     PetscInt ownership_range_lo;
00049     PetscInt ownership_range_hi;
00050     VecGetOwnershipRange(r_template, &ownership_range_lo, &ownership_range_hi);
00051     PetscInt local_size = ownership_range_hi - ownership_range_lo;
00052     PetscTools::SetupMat(mMassMatrix, 2*this->mpMesh->GetNumNodes(), 2*this->mpMesh->GetNumNodes(), 
00053                          2*this->mpMesh->CalculateMaximumNodeConnectivityPerProcess(),
00054                          local_size, local_size); 
00055 }
00056 
00057 
00058 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00059 void MatrixBasedBidomainSolver<ELEMENT_DIM,SPACE_DIM>::SetupLinearSystem(
00060         Vec currentSolution,
00061         bool computeMatrix)
00062 {
00063     assert(this->mpLinearSystem->rGetLhsMatrix() != NULL);
00064     assert(this->mpLinearSystem->rGetRhsVector() != NULL);
00065     assert(currentSolution != NULL);
00066 
00067     // create assembler
00068     if(!this->mpBidomainAssembler)
00069     {
00070         if(this->mBathSimulation)
00071         {
00072             this->mpBidomainAssembler = new BidomainWithBathAssembler<ELEMENT_DIM,SPACE_DIM>(this->mpMesh,this->mpBidomainTissue,this->mNumQuadPoints);
00073         }
00074         else
00075         {
00076             this->mpBidomainAssembler = new BidomainAssembler<ELEMENT_DIM,SPACE_DIM>(this->mpMesh,this->mpBidomainTissue,this->mNumQuadPoints);
00077         }
00078     }    
00079 
00080 
00082     // set up LHS matrix (and mass matrix)
00084     if (computeMatrix)
00085     {
00086         this->mpBidomainAssembler->SetMatrixToAssemble(this->mpLinearSystem->rGetLhsMatrix());
00087         this->mpBidomainAssembler->AssembleMatrix();
00088 
00089         // the BidomainMassMatrixAssembler deals with the mass matrix
00090         // for both bath and nonbath problems 
00091         assert(SPACE_DIM==ELEMENT_DIM);
00092         BidomainMassMatrixAssembler<SPACE_DIM> mass_matrix_assembler(this->mpMesh);
00093         mass_matrix_assembler.SetMatrixToAssemble(mMassMatrix);
00094         mass_matrix_assembler.Assemble();
00095 
00096         this->mpLinearSystem->AssembleIntermediateLhsMatrix();
00097         PetscMatTools::AssembleFinal(mMassMatrix);
00098     }
00099 
00100 
00101     HeartEventHandler::BeginEvent(HeartEventHandler::ASSEMBLE_RHS);
00102 
00104     // Set up z in b=Mz
00106     DistributedVectorFactory* p_factory = this->mpMesh->GetDistributedVectorFactory();
00107 
00108     // dist stripe for the current Voltage
00109     DistributedVector distributed_current_solution = p_factory->CreateDistributedVector(currentSolution);
00110     DistributedVector::Stripe distributed_current_solution_vm(distributed_current_solution, 0);
00111 
00112     // dist stripe for z
00113     DistributedVector dist_vec_matrix_based = p_factory->CreateDistributedVector(mVecForConstructingRhs);
00114     DistributedVector::Stripe dist_vec_matrix_based_vm(dist_vec_matrix_based, 0);
00115     DistributedVector::Stripe dist_vec_matrix_based_phie(dist_vec_matrix_based, 1);
00116 
00117     double Am = HeartConfig::Instance()->GetSurfaceAreaToVolumeRatio();
00118     double Cm  = HeartConfig::Instance()->GetCapacitance();
00119 
00120     if(!(this->mBathSimulation))
00121     {
00122         for (DistributedVector::Iterator index = dist_vec_matrix_based.Begin();
00123              index!= dist_vec_matrix_based.End();
00124              ++index)
00125         {
00126             double V = distributed_current_solution_vm[index];
00127             double F = - Am*this->mpBidomainTissue->rGetIionicCacheReplicated()[index.Global]
00128                        - this->mpBidomainTissue->rGetIntracellularStimulusCacheReplicated()[index.Global];
00129 
00130             dist_vec_matrix_based_vm[index] = Am*Cm*V*PdeSimulationTime::GetPdeTimeStepInverse() + F;
00131             dist_vec_matrix_based_phie[index] = 0.0;
00132         }
00133     }
00134     else
00135     {
00136         DistributedVector::Stripe dist_vec_matrix_based_phie(dist_vec_matrix_based, 1);
00137 
00138         for (DistributedVector::Iterator index = dist_vec_matrix_based.Begin();
00139              index!= dist_vec_matrix_based.End();
00140              ++index)
00141         {
00142     
00143             if ( !HeartRegionCode::IsRegionBath( this->mpMesh->GetNode(index.Global)->GetRegion() ))
00144             {
00145                 double V = distributed_current_solution_vm[index];
00146                 double F = - Am*this->mpBidomainTissue->rGetIionicCacheReplicated()[index.Global]
00147                            - this->mpBidomainTissue->rGetIntracellularStimulusCacheReplicated()[index.Global];
00148     
00149                 dist_vec_matrix_based_vm[index] = Am*Cm*V*PdeSimulationTime::GetPdeTimeStepInverse() + F;
00150             }
00151             else
00152             {
00153                 dist_vec_matrix_based_vm[index] = 0.0;
00154             }
00155     
00156             dist_vec_matrix_based_phie[index] = 0.0;
00157     
00158         }
00159     }
00160 
00161     dist_vec_matrix_based.Restore();
00162 
00164     // b = Mz
00166     MatMult(mMassMatrix, mVecForConstructingRhs, this->mpLinearSystem->rGetRhsVector());
00167 
00168     // assembling RHS is not finished yet, as Neumann bcs are added below, but
00169     // the event will be begun again inside this->mpBidomainAssembler->AssembleVector(); 
00170     HeartEventHandler::EndEvent(HeartEventHandler::ASSEMBLE_RHS);
00171 
00172 
00174     // apply Neumann boundary conditions
00176 
00177     this->mpBidomainAssembler->SetVectorToAssemble(this->mpLinearSystem->rGetRhsVector(), false/*don't zero vector!*/);
00178 
00179     // MUST be called every timestep in case the bcc has been reset, see comment in
00180     // this->ResetBoundaryConditionsContainer()
00181     this->mpBidomainAssembler->SetApplyNeummanBoundaryConditionsToVector(this->mpBoundaryConditions);
00182 
00183     this->mpBidomainAssembler->OnlyAssembleOnSurfaceElements();
00184     // note: don't need this for neumann bcs, would introduce parallel replication overhead
00185     //mpBidomainAssembler->SetCurrentSolution(currentSolution);
00186     this->mpBidomainAssembler->AssembleVector();
00187 
00188 
00190     // apply correction term
00192     if(mpBidomainCorrectionTermAssembler)
00193     {
00194         mpBidomainCorrectionTermAssembler->SetVectorToAssemble(this->mpLinearSystem->rGetRhsVector(), false/*don't zero vector!*/);
00195         // don't need to set current solution
00196         mpBidomainCorrectionTermAssembler->AssembleVector();
00197     }  
00198 
00199     this->mpLinearSystem->AssembleRhsVector();
00200 
00201     this->mpBoundaryConditions->ApplyDirichletToLinearProblem(*(this->mpLinearSystem), computeMatrix);
00202 
00203     if(this->mBathSimulation)
00204     {
00205         this->mpLinearSystem->AssembleFinalLhsMatrix();
00206         this->FinaliseForBath(computeMatrix,true);
00207     }
00208 
00209     if(computeMatrix)
00210     {
00211         this->mpLinearSystem->AssembleFinalLhsMatrix();
00212     }
00213     this->mpLinearSystem->AssembleRhsVector();
00214 }
00215 
00216 
00217 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00218 MatrixBasedBidomainSolver<ELEMENT_DIM,SPACE_DIM>::MatrixBasedBidomainSolver(
00219         bool bathSimulation,
00220         AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh,
00221         BidomainTissue<SPACE_DIM>* pTissue,
00222         BoundaryConditionsContainer<ELEMENT_DIM,SPACE_DIM,2>* pBoundaryConditions,
00223         unsigned numQuadPoints)
00224     : AbstractBidomainSolver<ELEMENT_DIM,SPACE_DIM>(bathSimulation,pMesh,pTissue,pBoundaryConditions)
00225 {
00226     // Tell tissue there's no need to replicate ionic caches
00227     pTissue->SetCacheReplication(false);
00228     mVecForConstructingRhs = NULL;
00229 
00230     if(HeartConfig::Instance()->GetUseStateVariableInterpolation())
00231     {
00232         mpBidomainCorrectionTermAssembler 
00233             = new BidomainCorrectionTermAssembler<ELEMENT_DIM,SPACE_DIM>(this->mpMesh,this->mpBidomainTissue,this->mNumQuadPoints);
00234         //We are going to need those caches after all
00235         pTissue->SetCacheReplication(true);
00236     }
00237     else
00238     {
00239         mpBidomainCorrectionTermAssembler = NULL;
00240     }
00241 }
00242 
00243 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM>
00244 MatrixBasedBidomainSolver<ELEMENT_DIM,SPACE_DIM>::~MatrixBasedBidomainSolver()
00245 {
00246     if(mVecForConstructingRhs)
00247     {
00248         VecDestroy(mVecForConstructingRhs);
00249         MatDestroy(mMassMatrix);
00250     }
00251     
00252     if(mpBidomainCorrectionTermAssembler)
00253     {
00254         delete mpBidomainCorrectionTermAssembler;
00255     }
00256 }
00257 
00259 // explicit instantiation
00261 
00262 template class MatrixBasedBidomainSolver<1,1>;
00263 template class MatrixBasedBidomainSolver<2,2>;
00264 template class MatrixBasedBidomainSolver<3,3>;
00265 

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