AbstractDynamicLinearPdeSolver.hpp

00001 
00002 /*
00003 
00004 Copyright (c) 2005-2015, University of Oxford.
00005 All rights reserved.
00006 
00007 University of Oxford means the Chancellor, Masters and Scholars of the
00008 University of Oxford, having an administrative office at Wellington
00009 Square, Oxford OX1 2JD, UK.
00010 
00011 This file is part of Chaste.
00012 
00013 Redistribution and use in source and binary forms, with or without
00014 modification, are permitted provided that the following conditions are met:
00015  * Redistributions of source code must retain the above copyright notice,
00016    this list of conditions and the following disclaimer.
00017  * Redistributions in binary form must reproduce the above copyright notice,
00018    this list of conditions and the following disclaimer in the documentation
00019    and/or other materials provided with the distribution.
00020  * Neither the name of the University of Oxford nor the names of its
00021    contributors may be used to endorse or promote products derived from this
00022    software without specific prior written permission.
00023 
00024 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
00025 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00026 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
00027 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
00028 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
00029 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
00030 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
00031 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
00032 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
00033 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00034 
00035 */
00036 
00037 #ifndef ABSTRACTDYNAMICLINEARPDESOLVER_HPP_
00038 #define ABSTRACTDYNAMICLINEARPDESOLVER_HPP_
00039 
00040 #include "Exception.hpp"
00041 #include "TimeStepper.hpp"
00042 #include "AbstractLinearPdeSolver.hpp"
00043 #include "PdeSimulationTime.hpp"
00044 #include "AbstractTimeAdaptivityController.hpp"
00045 #include "Hdf5DataWriter.hpp"
00046 #include "Hdf5ToVtkConverter.hpp"
00047 #include "Hdf5ToTxtConverter.hpp"
00048 
00055 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00056 class AbstractDynamicLinearPdeSolver : public AbstractLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>
00057 {
00058     friend class TestSimpleLinearParabolicSolver;
00059 protected:
00060 
00062     double mTstart;
00063 
00065     double mTend;
00066 
00068     bool mTimesSet;
00069 
00071     Vec mInitialCondition;
00072 
00074     bool mMatrixIsAssembled;
00075 
00080     bool mMatrixIsConstant;
00081 
00086     double mIdealTimeStep;
00087 
00089     double mLastWorkingTimeStep;
00090 
00092     AbstractTimeAdaptivityController* mpTimeAdaptivityController;
00093 
00098     bool mOutputToVtk;
00099 
00104     bool mOutputToParallelVtk;
00105 
00110     bool mOutputToTxt;
00111 
00113     std::string mOutputDirectory;
00114 
00116     std::string mFilenamePrefix;
00117 
00123     unsigned mPrintingTimestepMultiple;
00124 
00126     Hdf5DataWriter* mpHdf5Writer;
00127 
00129     std::vector<int> mVariableColumnIds;
00130 
00135     void InitialiseHdf5Writer();
00136 
00143     void WriteOneStep(double time, Vec solution);
00144 
00145 public:
00146 
00152     AbstractDynamicLinearPdeSolver(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh);
00153 
00160     void SetTimes(double tStart, double tEnd);
00161 
00167     void SetTimeStep(double dt);
00168 
00177     void SetInitialCondition(Vec initialCondition);
00178 
00182     virtual Vec Solve();
00183 
00185     void SetMatrixIsNotAssembled();
00186 
00192     void SetTimeAdaptivityController(AbstractTimeAdaptivityController* pTimeAdaptivityController);
00193 
00197     void SetOutputToVtk(bool output);
00198 
00202     void SetOutputToParallelVtk(bool output);
00203 
00207     void SetOutputToTxt(bool output);
00208 
00213     void SetOutputDirectoryAndPrefix(std::string outputDirectory, std::string prefix);
00214 
00219     void SetPrintingTimestepMultiple(unsigned multiple);
00220 };
00221 
00222 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00223 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::InitialiseHdf5Writer()
00224 {
00225     // Check that everything is set up correctly
00226     if ((mOutputDirectory=="") || (mFilenamePrefix==""))
00227     {
00228         EXCEPTION("Output directory or filename prefix has not been set");
00229     }
00230 
00231     // Create writer
00232     mpHdf5Writer = new Hdf5DataWriter(*(this->mpMesh)->GetDistributedVectorFactory(),
00233                                       mOutputDirectory,
00234                                       mFilenamePrefix);
00235 
00236     // Set writer to output all nodes
00237     mpHdf5Writer->DefineFixedDimension((this->mpMesh)->GetNumNodes());
00238 
00239     // Only used to get an estimate of the number of timesteps below
00240     unsigned estimated_num_printing_timesteps = 1u + (unsigned)((mTend - mTstart)/(mIdealTimeStep*mPrintingTimestepMultiple));
00241 
00248     assert(mVariableColumnIds.empty());
00249     for (unsigned i=0; i<PROBLEM_DIM; i++)
00250     {
00251         std::stringstream variable_name;
00252         variable_name << "Variable_" << i;
00253         mVariableColumnIds.push_back(mpHdf5Writer->DefineVariable(variable_name.str(),"undefined"));
00254     }
00255     mpHdf5Writer->DefineUnlimitedDimension("Time", "undefined", estimated_num_printing_timesteps);
00256 
00257     // End the define mode of the writer
00258     mpHdf5Writer->EndDefineMode();
00259 }
00260 
00261 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00262 AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::AbstractDynamicLinearPdeSolver(AbstractTetrahedralMesh<ELEMENT_DIM,SPACE_DIM>* pMesh)
00263     : AbstractLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>(pMesh),
00264       mTimesSet(false),
00265       mInitialCondition(NULL),
00266       mMatrixIsAssembled(false),
00267       mMatrixIsConstant(false),
00268       mIdealTimeStep(-1.0),
00269       mLastWorkingTimeStep(-1),
00270       mpTimeAdaptivityController(NULL),
00271       mOutputToVtk(false),
00272       mOutputToParallelVtk(false),
00273       mOutputToTxt(false),
00274       mOutputDirectory(""),
00275       mFilenamePrefix(""),
00276       mPrintingTimestepMultiple(1),
00277       mpHdf5Writer(NULL)
00278 {
00279 }
00280 
00281 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00282 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetTimes(double tStart, double tEnd)
00283 {
00284     mTstart = tStart;
00285     mTend = tEnd;
00286 
00287     if (mTstart >= mTend)
00288     {
00289         EXCEPTION("Start time has to be less than end time");
00290     }
00291 
00292     mTimesSet = true;
00293 }
00294 
00295 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00296 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetTimeStep(double dt)
00297 {
00298     if (dt <= 0)
00299     {
00300         EXCEPTION("Time step has to be greater than zero");
00301     }
00302 
00303     mIdealTimeStep = dt;
00304 }
00305 
00306 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00307 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetInitialCondition(Vec initialCondition)
00308 {
00309     assert(initialCondition != NULL);
00310     mInitialCondition = initialCondition;
00311 }
00312 
00313 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00314 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::WriteOneStep(double time, Vec solution)
00315 {
00316     mpHdf5Writer->PutUnlimitedVariable(time);
00317     if (PROBLEM_DIM == 1)
00318     {
00319         mpHdf5Writer->PutVector(mVariableColumnIds[0], solution);
00320     }
00321     else
00322     {
00323         mpHdf5Writer->PutStripedVector(mVariableColumnIds, solution);
00324     }
00325 }
00326 
00327 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00328 Vec AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::Solve()
00329 {
00330     // Begin by checking that everything has been set up correctly
00331     if (!mTimesSet)
00332     {
00333         EXCEPTION("SetTimes() has not been called");
00334     }
00335     if ((mIdealTimeStep <= 0.0) && (mpTimeAdaptivityController==NULL))
00336     {
00337         EXCEPTION("SetTimeStep() has not been called");
00338     }
00339     if (mInitialCondition == NULL)
00340     {
00341         EXCEPTION("SetInitialCondition() has not been called");
00342     }
00343 
00344     // If required, initialise HDF5 writer and output initial condition to HDF5 file
00345     bool print_output = (mOutputToVtk || mOutputToParallelVtk || mOutputToTxt);
00346     if (print_output)
00347     {
00348         InitialiseHdf5Writer();
00349         WriteOneStep(mTstart, mInitialCondition);
00350         mpHdf5Writer->AdvanceAlongUnlimitedDimension();
00351     }
00352 
00353     this->InitialiseForSolve(mInitialCondition);
00354 
00355     if (mIdealTimeStep < 0) // hasn't been set, so a controller must have been given
00356     {
00357         mIdealTimeStep = mpTimeAdaptivityController->GetNextTimeStep(mTstart, mInitialCondition);
00358     }
00359 
00360     /*
00361      * Note: we use the mIdealTimeStep here (the original timestep that was passed in, or
00362      * the last timestep suggested by the controller), rather than the last timestep used
00363      * (mLastWorkingTimeStep), because the timestep will be very slightly altered by the
00364      * stepper in the final timestep of the last printing-timestep-loop, and these floating
00365      * point errors can add up and eventually cause exceptions being thrown.
00366      */
00367     TimeStepper stepper(mTstart, mTend, mIdealTimeStep, mMatrixIsConstant);
00368 
00369     Vec solution = mInitialCondition;
00370     Vec next_solution;
00371 
00372     while (!stepper.IsTimeAtEnd())
00373     {
00374         bool timestep_changed = false;
00375 
00376         PdeSimulationTime::SetTime(stepper.GetTime());
00377 
00378         // Determine timestep to use
00379         double new_dt;
00380         if (mpTimeAdaptivityController)
00381         {
00382             // Get the timestep the controller wants to use and store it as the ideal timestep
00383             mIdealTimeStep = mpTimeAdaptivityController->GetNextTimeStep(stepper.GetTime(), solution);
00384 
00385             // Tell the stepper to use this timestep from now on...
00386             stepper.ResetTimeStep(mIdealTimeStep);
00387 
00388             // ..but now get the timestep from the stepper, as the stepper might need
00389             // to trim the timestep if it would take us over the end time
00390             new_dt = stepper.GetNextTimeStep();
00391             // Changes in timestep bigger than 0.001% will trigger matrix re-computation
00392             timestep_changed = (fabs(new_dt/mLastWorkingTimeStep - 1.0) > 1e-5);
00393         }
00394         else
00395         {
00396             new_dt = stepper.GetNextTimeStep();
00397 
00398             //new_dt should be roughly the same size as mIdealTimeStep - we should never need to take a tiny step
00399 
00400             if (mMatrixIsConstant && fabs(new_dt/mIdealTimeStep - 1.0) > 1e-5)
00401             {
00402                 // Here we allow for changes of up to 0.001%
00403                 // Note that the TimeStepper guarantees that changes in dt are no bigger than DBL_EPSILON*current_time
00404                 NEVER_REACHED;
00405             }
00406         }
00407 
00408         // Save the timestep as the last one use, and also put it in PdeSimulationTime
00409         // so everyone can see it
00410         mLastWorkingTimeStep = new_dt;
00411         PdeSimulationTime::SetPdeTimeStepAndNextTime(new_dt, stepper.GetNextTime());
00412 
00413         // Solve
00414         try
00415         {
00416             // (This runs the cell ODE models in heart simulations)
00417             this->PrepareForSetupLinearSystem(solution);
00418         }
00419         catch(Exception& e)
00420         {
00421             // We only need to clean up memory if we are NOT on the first PDE time step,
00422             // as someone else cleans up the mInitialCondition vector in higher classes.
00423             if (solution != mInitialCondition)
00424             {
00425                 HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
00426                 PetscTools::Destroy(solution);
00427                 HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
00428             }
00429             throw e;
00430         }
00431 
00432         bool compute_matrix = (!mMatrixIsConstant || !mMatrixIsAssembled || timestep_changed);
00433 
00434         this->SetupLinearSystem(solution, compute_matrix);
00435 
00436         this->FinaliseLinearSystem(solution);
00437 
00438         if (compute_matrix)
00439         {
00440             this->mpLinearSystem->ResetKspSolver();
00441         }
00442 
00443         next_solution = this->mpLinearSystem->Solve(solution);
00444 
00445         if (mMatrixIsConstant)
00446         {
00447             mMatrixIsAssembled = true;
00448         }
00449 
00450         this->FollowingSolveLinearSystem(next_solution);
00451 
00452         stepper.AdvanceOneTimeStep();
00453 
00454         // Avoid memory leaks
00455         if (solution != mInitialCondition)
00456         {
00457             HeartEventHandler::BeginEvent(HeartEventHandler::COMMUNICATION);
00458             PetscTools::Destroy(solution);
00459             HeartEventHandler::EndEvent(HeartEventHandler::COMMUNICATION);
00460         }
00461         solution = next_solution;
00462 
00463         // If required, output next solution to HDF5 file
00464         if (print_output && (stepper.GetTotalTimeStepsTaken()%mPrintingTimestepMultiple == 0) )
00465         {
00466             WriteOneStep(stepper.GetTime(), solution);
00467             mpHdf5Writer->AdvanceAlongUnlimitedDimension();
00468         }
00469     }
00470 
00471     // Avoid memory leaks
00472     if (mpHdf5Writer != NULL)
00473     {
00474         delete mpHdf5Writer;
00475         mpHdf5Writer = NULL;
00476     }
00477 
00478     // Convert HDF5 output to other formats as required
00479 
00480     if (mOutputToVtk)
00481     {
00482         Hdf5ToVtkConverter<ELEMENT_DIM,SPACE_DIM> converter(FileFinder(mOutputDirectory, RelativeTo::ChasteTestOutput),
00483                                                             mFilenamePrefix, this->mpMesh, false, false);
00484     }
00485     if (mOutputToParallelVtk)
00486     {
00487         Hdf5ToVtkConverter<ELEMENT_DIM,SPACE_DIM> converter(FileFinder(mOutputDirectory, RelativeTo::ChasteTestOutput),
00488                                                             mFilenamePrefix, this->mpMesh, true, false);
00489     }
00490     if (mOutputToTxt)
00491     {
00492         Hdf5ToTxtConverter<ELEMENT_DIM,SPACE_DIM> converter(FileFinder(mOutputDirectory, RelativeTo::ChasteTestOutput),
00493                                                             mFilenamePrefix, this->mpMesh);
00494     }
00495 
00496     return solution;
00497 }
00498 
00499 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00500 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetMatrixIsNotAssembled()
00501 {
00502     mMatrixIsAssembled = false;
00503 }
00504 
00505 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00506 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetTimeAdaptivityController(AbstractTimeAdaptivityController* pTimeAdaptivityController)
00507 {
00508     assert(pTimeAdaptivityController != NULL);
00509     assert(mpTimeAdaptivityController == NULL);
00510     mpTimeAdaptivityController = pTimeAdaptivityController;
00511 }
00512 
00513 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00514 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetOutputToVtk(bool output)
00515 {
00516     mOutputToVtk = output;
00517 }
00518 
00519 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00520 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetOutputToParallelVtk(bool output)
00521 {
00522     mOutputToParallelVtk = output;
00523 }
00524 
00525 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00526 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetOutputToTxt(bool output)
00527 {
00528     mOutputToTxt = output;
00529 }
00530 
00531 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00532 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetOutputDirectoryAndPrefix(std::string outputDirectory, std::string prefix)
00533 {
00534     mOutputDirectory = outputDirectory;
00535     mFilenamePrefix = prefix;
00536 }
00537 
00538 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM>
00539 void AbstractDynamicLinearPdeSolver<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM>::SetPrintingTimestepMultiple(unsigned multiple)
00540 {
00541     mPrintingTimestepMultiple = multiple;
00542 }
00543 
00544 #endif /*ABSTRACTDYNAMICLINEARPDESOLVER_HPP_*/

Generated by  doxygen 1.6.2