AbstractDynamicLinearPdeSolver.hpp

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 #ifndef ABSTRACTDYNAMICLINEARPDESOLVER_HPP_
00031 #define ABSTRACTDYNAMICLINEARPDESOLVER_HPP_
00032 
00033 #include "TimeStepper.hpp"
00034 #include "AbstractLinearPdeSolver.hpp"
00035 #include "PdeSimulationTime.hpp"
00036 #include "AbstractTimeAdaptivityController.hpp"
00037 
00038 
00045 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00046 class AbstractDynamicLinearPdeSolver : public AbstractLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>
00047 {
00048     friend class TestSimpleLinearParabolicSolver;
00049 
00050 protected:
00052     double mTstart;
00053 
00055     double mTend;
00056 
00058     bool mTimesSet;
00059 
00061     Vec mInitialCondition;
00062 
00064     bool mMatrixIsAssembled;
00065 
00068     bool mMatrixIsConstant;
00069 
00073     double mIdealTimeStep;
00074     
00076     double mLastWorkingTimeStep;
00077 
00079     AbstractTimeAdaptivityController* mpTimeAdaptivityController;
00080 
00081 public:
00086     AbstractDynamicLinearPdeSolver(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh);
00087 
00093     void SetTimes(double tStart, double tEnd);
00094 
00098     void SetTimeStep(double dt);
00099 
00104     void SetInitialCondition(Vec initialCondition);
00105 
00107     Vec Solve();
00108 
00110     void SetMatrixIsNotAssembled()
00111     {
00112         mMatrixIsAssembled = false;
00113     }
00114     
00118     void SetTimeAdaptivityController(AbstractTimeAdaptivityController* pTimeAdaptivityController)
00119     {
00120         assert(pTimeAdaptivityController != NULL);
00121         assert(mpTimeAdaptivityController == NULL);
00122         mpTimeAdaptivityController = pTimeAdaptivityController;
00123     }
00124 };    
00125 
00126 
00127 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00128 AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::AbstractDynamicLinearPdeSolver(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh)
00129     : AbstractLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>(pMesh),
00130       mTimesSet(false),
00131       mInitialCondition(NULL),
00132       mMatrixIsAssembled(false),
00133       mMatrixIsConstant(false),
00134       mIdealTimeStep(-1.0),
00135       mLastWorkingTimeStep(-1),
00136       mpTimeAdaptivityController(NULL)
00137 {
00138 }
00139 
00140 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00141 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetTimes(double tStart, double tEnd)
00142 {
00143     mTstart = tStart;
00144     mTend   = tEnd;
00145 
00146     if (mTstart >= mTend)
00147     {
00148         EXCEPTION("Starting time has to less than ending time");
00149     }
00150 
00151     mTimesSet = true;
00152 }
00153 
00154 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00155 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetTimeStep(double dt)
00156 {
00157     if (dt <= 0)
00158     {
00159         EXCEPTION("Time step has to be greater than zero");
00160     }
00161 
00162     mIdealTimeStep = dt;
00163 }
00164 
00165 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00166 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetInitialCondition(Vec initialCondition)
00167 {
00168     assert(initialCondition!=NULL);
00169     mInitialCondition = initialCondition;
00170 }
00171 
00172 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00173 Vec AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::Solve()
00174 {
00175     assert(mTimesSet);
00176     assert(mIdealTimeStep > 0.0 || mpTimeAdaptivityController);
00177     assert(mInitialCondition != NULL);
00178 
00179     this->InitialiseForSolve(mInitialCondition);
00180 
00181     if(mIdealTimeStep < 0) // hasn't been set, so a controller must have been given
00182     {
00183         mIdealTimeStep = mpTimeAdaptivityController->GetNextTimeStep(mTstart, mInitialCondition);
00184     }
00185 
00186     // Note: we use the mIdealTimeStep here (the original timestep that was passed in, or
00187     // the last timestep suggested by the controller), rather than the last timestep used
00188     // (mLastWorkingTimeStep), because the timestep will be very slightly altered by
00189     // the stepper in the final timestep of the last printing-timestep-loop, and these
00190     // floating point errors can add up and eventually cause exceptions being thrown.
00191     TimeStepper stepper(mTstart, mTend, mIdealTimeStep, mMatrixIsConstant);
00192 
00193     Vec current_solution = mInitialCondition;
00194     Vec next_solution;
00195 
00196     while ( !stepper.IsTimeAtEnd() )
00197     {    
00198         bool timestep_changed = false;
00199 
00200         PdeSimulationTime::SetTime(stepper.GetTime());
00201         
00203         // determine timestep to use 
00205         double new_dt;
00206         if (mpTimeAdaptivityController)
00207         {
00208             // get the timestep the controller wants to use and store it as the ideal timestep
00209             mIdealTimeStep = mpTimeAdaptivityController->GetNextTimeStep(stepper.GetTime(), current_solution);
00210             // tell the stepper to use this timestep from now on.
00211             stepper.ResetTimeStep(mIdealTimeStep);
00212             // ..but now get the timestep from the stepper, as the stepper might need
00213             // to trim the timestep if it would take us over the end time
00214             new_dt = stepper.GetNextTimeStep();
00215 
00216             timestep_changed = (fabs(mLastWorkingTimeStep-new_dt) > 1e-8);
00217         }
00218         else
00219         {
00220             new_dt = stepper.GetNextTimeStep();
00221         }
00222 
00223         // save the timestep as the last one use, and also put it in PdeSimulationTime
00224         // so everyone can see it
00225         mLastWorkingTimeStep = new_dt;
00226         PdeSimulationTime::SetPdeTimeStep( new_dt ); 
00227 
00229         // solve 
00231 
00232         this->PrepareForSetupLinearSystem(current_solution);
00233 
00234         bool compute_matrix = (!mMatrixIsConstant || !mMatrixIsAssembled || timestep_changed);
00235 //        if (compute_matrix) std::cout << " ** ASSEMBLING MATRIX!!! ** " << std::endl;
00236         this->SetupLinearSystem(current_solution, compute_matrix);
00237        
00238         this->FinaliseLinearSystem(current_solution);
00239         
00240         if (compute_matrix)
00241         {
00242             this->mpLinearSystem->ResetKspSolver();
00243         }
00244     
00245         next_solution = this->mpLinearSystem->Solve(current_solution);
00246 
00247         if (mMatrixIsConstant)
00248         {
00249             mMatrixIsAssembled = true;
00250         }
00251 
00252         this->FollowingSolveLinearSystem(next_solution);
00253 
00254         stepper.AdvanceOneTimeStep();
00255 
00256         // Avoid memory leaks
00257         if (current_solution != mInitialCondition)
00258         {
00259             HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
00260             VecDestroy(current_solution);
00261             HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
00262         }
00263         current_solution = next_solution;
00264     }    
00265     
00266     return current_solution;
00267 }
00268 
00269 
00270 #endif /*ABSTRACTDYNAMICLINEARPDESOLVER_HPP_*/

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