00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
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
00226 if ((mOutputDirectory=="") || (mFilenamePrefix==""))
00227 {
00228 EXCEPTION("Output directory or filename prefix has not been set");
00229 }
00230
00231
00232 mpHdf5Writer = new Hdf5DataWriter(*(this->mpMesh)->GetDistributedVectorFactory(),
00233 mOutputDirectory,
00234 mFilenamePrefix);
00235
00236
00237 mpHdf5Writer->DefineFixedDimension((this->mpMesh)->GetNumNodes());
00238
00239
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
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
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
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)
00356 {
00357 mIdealTimeStep = mpTimeAdaptivityController->GetNextTimeStep(mTstart, mInitialCondition);
00358 }
00359
00360
00361
00362
00363
00364
00365
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
00379 double new_dt;
00380 if (mpTimeAdaptivityController)
00381 {
00382
00383 mIdealTimeStep = mpTimeAdaptivityController->GetNextTimeStep(stepper.GetTime(), solution);
00384
00385
00386 stepper.ResetTimeStep(mIdealTimeStep);
00387
00388
00389
00390 new_dt = stepper.GetNextTimeStep();
00391
00392 timestep_changed = (fabs(new_dt/mLastWorkingTimeStep - 1.0) > 1e-5);
00393 }
00394 else
00395 {
00396 new_dt = stepper.GetNextTimeStep();
00397
00398
00399
00400 if (mMatrixIsConstant && fabs(new_dt/mIdealTimeStep - 1.0) > 1e-5)
00401 {
00402
00403
00404 NEVER_REACHED;
00405 }
00406 }
00407
00408
00409
00410 mLastWorkingTimeStep = new_dt;
00411 PdeSimulationTime::SetPdeTimeStepAndNextTime(new_dt, stepper.GetNextTime());
00412
00413
00414 try
00415 {
00416
00417 this->PrepareForSetupLinearSystem(solution);
00418 }
00419 catch(Exception& e)
00420 {
00421
00422
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
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
00464 if (print_output && (stepper.GetTotalTimeStepsTaken()%mPrintingTimestepMultiple == 0) )
00465 {
00466 WriteOneStep(stepper.GetTime(), solution);
00467 mpHdf5Writer->AdvanceAlongUnlimitedDimension();
00468 }
00469 }
00470
00471
00472 if (mpHdf5Writer != NULL)
00473 {
00474 delete mpHdf5Writer;
00475 mpHdf5Writer = NULL;
00476 }
00477
00478
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