35#ifndef LINEARPARABOLICPDESYSTEMWITHCOUPLEDODESYSTEMSOLVER_HPP_
36#define LINEARPARABOLICPDESYSTEMWITHCOUPLEDODESYSTEMSOLVER_HPP_
38#include "AbstractAssemblerSolverHybrid.hpp"
39#include "AbstractDynamicLinearPdeSolver.hpp"
40#include "AbstractLinearParabolicPdeSystemForCoupledOdeSystem.hpp"
41#include "TetrahedralMesh.hpp"
42#include "BoundaryConditionsContainer.hpp"
43#include "AbstractOdeSystemForCoupledPdeSystem.hpp"
44#include "CvodeAdaptor.hpp"
45#include "BackwardEulerIvpOdeSolver.hpp"
46#include "Warnings.hpp"
47#include "VtkMeshWriter.hpp"
49#include <boost/shared_ptr.hpp>
60template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM=ELEMENT_DIM,
unsigned PROBLEM_DIM=1>
117 c_vector<double, ELEMENT_DIM+1>& rPhi,
118 c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
120 c_vector<double,PROBLEM_DIM>& rU,
121 c_matrix<double, PROBLEM_DIM, SPACE_DIM>& rGradU,
135 c_vector<double, ELEMENT_DIM+1>& rPhi,
136 c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
138 c_vector<double,PROBLEM_DIM>& rU,
139 c_matrix<double,PROBLEM_DIM,SPACE_DIM>& rGradU,
190 std::vector<AbstractOdeSystemForCoupledPdeSystem*> odeSystemsAtNodes=std::vector<AbstractOdeSystemForCoupledPdeSystem*>(),
191 boost::shared_ptr<AbstractIvpOdeSolver> pOdeSolver=boost::shared_ptr<AbstractIvpOdeSolver>());
250template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
252 c_vector<double, ELEMENT_DIM+1>& rPhi,
253 c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
255 c_vector<double,PROBLEM_DIM>& rU,
256 c_matrix<double, PROBLEM_DIM, SPACE_DIM>& rGradU,
260 c_matrix<
double, PROBLEM_DIM*(ELEMENT_DIM+1), PROBLEM_DIM*(ELEMENT_DIM+1)> matrix_term = zero_matrix<double>(PROBLEM_DIM*(ELEMENT_DIM+1), PROBLEM_DIM*(ELEMENT_DIM+1));
263 for (
unsigned pde_index=0; pde_index<PROBLEM_DIM; pde_index++)
266 c_matrix<double, SPACE_DIM, SPACE_DIM> this_pde_diffusion_term = mpPdeSystem->
ComputeDiffusionTerm(rX, pde_index, pElement);
267 c_matrix<
double, 1*(ELEMENT_DIM+1), 1*(ELEMENT_DIM+1)> this_stiffness_matrix =
268 prod(trans(rGradPhi), c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>(prod(this_pde_diffusion_term, rGradPhi)) )
269 + timestep_inverse * this_dudt_coefficient * outer_prod(rPhi, rPhi);
271 for (
unsigned i=0; i<ELEMENT_DIM+1; i++)
273 for (
unsigned j=0; j<ELEMENT_DIM+1; j++)
275 matrix_term(i*PROBLEM_DIM + pde_index, j*PROBLEM_DIM + pde_index) = this_stiffness_matrix(i,j);
282template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
284 c_vector<double, ELEMENT_DIM+1>& rPhi,
285 c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
287 c_vector<double,PROBLEM_DIM>& rU,
288 c_matrix<double,PROBLEM_DIM,SPACE_DIM>& rGradU,
292 c_vector<
double, PROBLEM_DIM*(ELEMENT_DIM+1)> vector_term;
293 vector_term = zero_vector<double>(PROBLEM_DIM*(ELEMENT_DIM+1));
296 for (
unsigned pde_index=0; pde_index<PROBLEM_DIM; pde_index++)
299 double this_source_term = mpPdeSystem->
ComputeSourceTerm(rX, rU, mInterpolatedOdeStateVariables, pde_index);
300 c_vector<double, ELEMENT_DIM+1> this_vector_term;
301 this_vector_term = (this_source_term + timestep_inverse*this_dudt_coefficient*rU(pde_index))* rPhi;
303 for (
unsigned i=0; i<ELEMENT_DIM+1; i++)
305 vector_term(i*PROBLEM_DIM + pde_index) = this_vector_term(i);
312template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
315 mInterpolatedOdeStateVariables.clear();
317 if (mOdeSystemsPresent)
319 unsigned num_state_variables = mOdeSystemsAtNodes[0]->GetNumberOfStateVariables();
320 mInterpolatedOdeStateVariables.resize(num_state_variables, 0.0);
324template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
327 if (mOdeSystemsPresent)
329 unsigned num_state_variables = mOdeSystemsAtNodes[0]->GetNumberOfStateVariables();
331 for (
unsigned i=0; i<num_state_variables; i++)
333 mInterpolatedOdeStateVariables[i] += phiI * mOdeSystemsAtNodes[pNode->
GetIndex()]->rGetStateVariables()[i];
338template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
341 if (this->mpLinearSystem == NULL)
349 preallocation *= PROBLEM_DIM;
356 this->mpLinearSystem =
new LinearSystem(initialSolution, preallocation);
359 assert(this->mpLinearSystem);
360 this->mpLinearSystem->SetMatrixIsSymmetric(
true);
361 this->mpLinearSystem->SetKspType(
"cg");
364template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
367 this->SetupGivenLinearSystem(currentSolution, computeMatrix, this->mpLinearSystem);
370template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
375 std::vector<AbstractOdeSystemForCoupledPdeSystem*> odeSystemsAtNodes,
376 boost::shared_ptr<AbstractIvpOdeSolver> pOdeSolver)
380 mpPdeSystem(pPdeSystem),
381 mOdeSystemsAtNodes(odeSystemsAtNodes),
382 mpOdeSolver(pOdeSolver),
384 mOdeSystemsPresent(false),
385 mClearOutputDirectory(false)
414template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
417 if (mOdeSystemsPresent)
419 for (
unsigned i=0; i<mOdeSystemsAtNodes.size(); i++)
421 delete mOdeSystemsAtNodes[i];
426template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
429 if (mOdeSystemsPresent)
436 std::vector<double> current_soln_this_node(PROBLEM_DIM);
439 for (
unsigned node_index=0; node_index<mpMesh->
GetNumNodes(); node_index++)
442 for (
unsigned pde_index=0; pde_index<PROBLEM_DIM; pde_index++)
444 double current_soln_this_pde_this_node = soln_repl[PROBLEM_DIM*node_index + pde_index];
445 current_soln_this_node[pde_index] = current_soln_this_pde_this_node;
449 mOdeSystemsAtNodes[node_index]->SetPdeSolution(current_soln_this_node);
452 mpOdeSolver->SolveAndUpdateStateVariable(mOdeSystemsAtNodes[node_index], time, next_time, dt);
457template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
460 mClearOutputDirectory = clearDirectory;
461 this->mOutputDirectory = outputDirectory;
464template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
467 assert(samplingTimeStep >= this->mIdealTimeStep);
468 mSamplingTimeStep = samplingTimeStep;
471template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
475 if (this->mOutputDirectory ==
"")
477 EXCEPTION(
"SetOutputDirectory() must be called prior to SolveAndWriteResultsToFile()");
479 if (this->mTimesSet ==
false)
481 EXCEPTION(
"SetTimes() must be called prior to SolveAndWriteResultsToFile()");
483 if (this->mIdealTimeStep <= 0.0)
485 EXCEPTION(
"SetTimeStep() must be called prior to SolveAndWriteResultsToFile()");
489 EXCEPTION(
"SetSamplingTimeStep() must be called prior to SolveAndWriteResultsToFile()");
491 if (!this->mInitialCondition)
493 EXCEPTION(
"SetInitialCondition() must be called prior to SolveAndWriteResultsToFile()");
498 OutputFileHandler output_file_handler(this->mOutputDirectory, mClearOutputDirectory);
499 mpVtkMetaFile = output_file_handler.
OpenOutputFile(
"results.pvd");
500 *mpVtkMetaFile <<
"<?xml version=\"1.0\"?>\n";
501 *mpVtkMetaFile <<
"<VTKFile type=\"Collection\" version=\"0.1\" byte_order=\"LittleEndian\" compressor=\"vtkZLibDataCompressor\">\n";
502 *mpVtkMetaFile <<
" <Collection>\n";
505 Vec initial_condition = this->mInitialCondition;
506 WriteVtkResultsToFile(initial_condition, 0);
509 TimeStepper stepper(this->mTstart, this->mTend, mSamplingTimeStep);
518 Vec soln = this->Solve();
521 if (this->mInitialCondition != initial_condition)
525 this->mInitialCondition = soln;
535 if (this->mInitialCondition != initial_condition)
539 this->mInitialCondition = initial_condition;
542 *mpVtkMetaFile <<
" </Collection>\n";
543 *mpVtkMetaFile <<
"</VTKFile>\n";
544 mpVtkMetaFile->close();
547 WARNING(
"VTK is not installed and is required for this functionality");
552template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
558 std::stringstream time;
559 time << numTimeStepsElapsed;
569 for (
unsigned pde_index=0; pde_index<PROBLEM_DIM; pde_index++)
572 std::vector<double> pde_index_data;
573 pde_index_data.resize(num_nodes, 0.0);
574 for (
unsigned node_index=0; node_index<num_nodes; node_index++)
576 pde_index_data[node_index] = solution_repl[PROBLEM_DIM*node_index + pde_index];
580 std::stringstream data_name;
581 data_name <<
"PDE variable " << pde_index;
582 mesh_writer.
AddPointData(data_name.str(), pde_index_data);
585 if (mOdeSystemsPresent)
593 std::vector<std::vector<double> > ode_data;
594 unsigned num_odes = mOdeSystemsAtNodes[0]->rGetStateVariables().size();
595 ode_data.resize(num_odes);
596 for (
unsigned ode_index=0; ode_index<num_odes; ode_index++)
598 ode_data[ode_index].resize(num_nodes, 0.0);
601 for (
unsigned node_index=0; node_index<num_nodes; node_index++)
603 std::vector<double> all_odes_this_node = mOdeSystemsAtNodes[node_index]->rGetStateVariables();
604 for (
unsigned i=0; i<num_odes; i++)
606 ode_data[i][node_index] = all_odes_this_node[i];
610 for (
unsigned ode_index=0; ode_index<num_odes; ode_index++)
612 std::vector<double> ode_index_data = ode_data[ode_index];
615 std::stringstream data_name;
616 data_name <<
"ODE variable " << ode_index;
617 mesh_writer.
AddPointData(data_name.str(), ode_index_data);
622 *mpVtkMetaFile <<
" <DataSet timestep=\"";
623 *mpVtkMetaFile << numTimeStepsElapsed;
624 *mpVtkMetaFile <<
"\" group=\"\" part=\"0\" file=\"results_";
625 *mpVtkMetaFile << numTimeStepsElapsed;
626 *mpVtkMetaFile <<
".vtu\"/>\n";
630template<
unsigned ELEMENT_DIM,
unsigned SPACE_DIM,
unsigned PROBLEM_DIM>
633 return mOdeSystemsAtNodes[index];
const double DOUBLE_UNSET
#define EXCEPTION(message)
BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > * mpBoundaryConditions
virtual double ComputeSourceTerm(const ChastePoint< SPACE_DIM > &rX, c_vector< double, PROBLEM_DIM > &rU, std::vector< double > &rOdeSolution, unsigned pdeIndex)=0
virtual c_matrix< double, SPACE_DIM, SPACE_DIM > ComputeDiffusionTerm(const ChastePoint< SPACE_DIM > &rX, unsigned pdeIndex, Element< ELEMENT_DIM, SPACE_DIM > *pElement=NULL)=0
virtual double ComputeDuDtCoefficientFunction(const ChastePoint< SPACE_DIM > &rX, unsigned pdeIndex)=0
virtual unsigned GetNumNodes() const
unsigned CalculateMaximumContainingElementsPerProcess() const
void SolveAndWriteResultsToFile()
AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > * mpMesh
LinearParabolicPdeSystemWithCoupledOdeSystemSolver(TetrahedralMesh< ELEMENT_DIM, SPACE_DIM > *pMesh, AbstractLinearParabolicPdeSystemForCoupledOdeSystem< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > *pPdeSystem, BoundaryConditionsContainer< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > *pBoundaryConditions, std::vector< AbstractOdeSystemForCoupledPdeSystem * > odeSystemsAtNodes=std::vector< AbstractOdeSystemForCoupledPdeSystem * >(), boost::shared_ptr< AbstractIvpOdeSolver > pOdeSolver=boost::shared_ptr< AbstractIvpOdeSolver >())
void ResetInterpolatedQuantities()
void IncrementInterpolatedQuantities(double phiI, const Node< SPACE_DIM > *pNode)
std::vector< AbstractOdeSystemForCoupledPdeSystem * > mOdeSystemsAtNodes
void PrepareForSetupLinearSystem(Vec currentPdeSolution)
boost::shared_ptr< AbstractIvpOdeSolver > mpOdeSolver
void SetOutputDirectory(std::string outputDirectory, bool clearDirectory=false)
void SetupLinearSystem(Vec currentSolution, bool computeMatrix)
c_matrix< double, PROBLEM_DIM *(ELEMENT_DIM+1), PROBLEM_DIM *(ELEMENT_DIM+1)> ComputeMatrixTerm(c_vector< double, ELEMENT_DIM+1 > &rPhi, c_matrix< double, SPACE_DIM, ELEMENT_DIM+1 > &rGradPhi, ChastePoint< SPACE_DIM > &rX, c_vector< double, PROBLEM_DIM > &rU, c_matrix< double, PROBLEM_DIM, SPACE_DIM > &rGradU, Element< ELEMENT_DIM, SPACE_DIM > *pElement)
c_vector< double, PROBLEM_DIM *(ELEMENT_DIM+1)> ComputeVectorTerm(c_vector< double, ELEMENT_DIM+1 > &rPhi, c_matrix< double, SPACE_DIM, ELEMENT_DIM+1 > &rGradPhi, ChastePoint< SPACE_DIM > &rX, c_vector< double, PROBLEM_DIM > &rU, c_matrix< double, PROBLEM_DIM, SPACE_DIM > &rGradU, Element< ELEMENT_DIM, SPACE_DIM > *pElement)
void SetSamplingTimeStep(double samplingTimeStep)
std::vector< double > mInterpolatedOdeStateVariables
AbstractOdeSystemForCoupledPdeSystem * GetOdeSystemAtNode(unsigned index)
void InitialiseForSolve(Vec initialSolution=NULL)
~LinearParabolicPdeSystemWithCoupledOdeSystemSolver()
void WriteVtkResultsToFile()
AbstractLinearParabolicPdeSystemForCoupledOdeSystem< ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM > * mpPdeSystem
bool mClearOutputDirectory
unsigned GetIndex() const
out_stream OpenOutputFile(const std::string &rFileName, std::ios_base::openmode mode=std::ios::out|std::ios::trunc) const
static double GetPdeTimeStep()
static double GetPdeTimeStepInverse()
static double GetNextTime()
unsigned GetTotalTimeStepsTaken() const
void AdvanceOneTimeStep()
double GetNextTime() const
void AddPointData(std::string name, std::vector< double > data)
void WriteFilesUsingMesh(AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM > &rMesh, bool keepOriginalElementIndexing=true)