#ifndef eulerPROBLEM_IMPL_H_
#define eulerPROBLEM_IMPL_H_
#include <core/mfilename.h>
#include <matrices/tnlMatrixSetter.h>
#include <solvers/pde/tnlExplicitUpdater.h>
#include <solvers/pde/tnlLinearSystemAssembler.h>
#include <solvers/pde/tnlBackwardTimeDiscretisation.h>
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
tnlString
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
getTypeStatic()
{
return tnlString( "eulerProblem< " ) + Mesh :: getTypeStatic() + " >";
}
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
tnlString
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
getPrologHeader() const
{
return tnlString( "euler" );
}
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
void
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
writeProlog( tnlLogger& logger, const tnlParameterContainer& parameters ) const
{
/****
* Add data you want to have in the computation report (log) as follows:
* logger.writeParameter< double >( "Parameter description", parameter );
*/
}
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
bool
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
setup( const tnlParameterContainer& parameters )
{
if( ! this->boundaryCondition.setup( parameters, "boundary-conditions-" ) ||
! this->rightHandSide.setup( parameters, "right-hand-side-" ) )
return false;
return true;
}
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
typename eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::IndexType
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
getDofs( const MeshType& mesh ) const
{
/****
* Return number of DOFs (degrees of freedom) i.e. number
* of unknowns to be resolved by the main solver.
*/
return 3*mesh.template getEntitiesCount< typename MeshType::Cell >();
}
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
void
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
bindDofs( const MeshType& mesh,
DofVectorType& dofVector )
{
}
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
bool
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
setInitialCondition( const tnlParameterContainer& parameters,
const MeshType& mesh,
DofVectorType& dofs,
MeshDependentDataType& meshDependentData )
{
typedef typename MeshType::Cell Cell;
double gamma = parameters.getParameter< double >( "gamma" );
double rhoL = parameters.getParameter< double >( "left-density" );
double velL = parameters.getParameter< double >( "left-velocity" );
double preL = parameters.getParameter< double >( "left-pressure" );
double eL = ( preL / (gamma - 1) ) + 0.5 * rhoL * velL * velL;
double rhoR = parameters.getParameter< double >( "right-density" );
double velR = parameters.getParameter< double >( "right-velocity" );
double preR = parameters.getParameter< double >( "right-pressure" );
double eR = ( preR / (gamma - 1) ) + 0.5 * rhoR * velR * velR;
double x0 = parameters.getParameter< double >( "riemann-border" );
cout << gamma << " " << rhoL << " " << velL << " " << preL << " " << eL << " " << rhoR << " " << velR << " " << preR << " " << eR << " " << x0 << " " << gamma << endl;
int count = mesh.template getEntitiesCount< Cell >();
this->rho.bind(dofs,0,count);
this->rhoVel.bind(dofs,count,count);
this->energy.bind(dofs,2 * count,count);
this->data.setSize(2*count);
this->pressure.bind(this->data,0,count);
this->velocity.bind(this->data,count,count);
for(long int i = 0; i < count; i++)
if (i < x0 * count )
{
this->rho[i] = rhoL;
this->rhoVel[i] = rhoL * velL;
this->energy[i] = eL;
this->velocity[i] = velL;
this->pressure[i] = preL;
}
else
{
this->rho[i] = rhoR;
this->rhoVel[i] = rhoR * velR;
this->energy[i] = eR;
this->velocity[i] = velR;
this->pressure[i] = preR;
};
this->gamma = gamma;
cout << "dofs = " << dofs << endl;
getchar();
/*
const tnlString& initialConditionFile = parameters.getParameter< tnlString >( "initial-condition" );
if( ! dofs.load( initialConditionFile ) )
{
cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." << endl;
return false;
}
*/
return true;
}
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
template< typename Matrix >
bool
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
setupLinearSystem( const MeshType& mesh,
Matrix& matrix )
{
const IndexType dofs = this->getDofs( mesh );
typedef typename Matrix::CompressedRowsLengthsVector CompressedRowsLengthsVectorType;
CompressedRowsLengthsVectorType rowLengths;
if( ! rowLengths.setSize( dofs ) )
return false;
tnlMatrixSetter< MeshType, DifferentialOperator, BoundaryCondition, CompressedRowsLengthsVectorType > matrixSetter;
matrixSetter.template getCompressedRowsLengths< typename Mesh::Cell >( mesh,
differentialOperator,
boundaryCondition,
rowLengths );
matrix.setDimensions( dofs, dofs );
if( ! matrix.setCompressedRowsLengths( rowLengths ) )
return false;
return true;
}
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
bool
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
makeSnapshot( const RealType& time,
const IndexType& step,
const MeshType& mesh,
DofVectorType& dofs,
MeshDependentDataType& meshDependentData )
{
cout << endl << "Writing output at time " << time << " step " << step << "." << endl;
this->bindDofs( mesh, dofs );
tnlString fileName;
ofstream vysledek;
cout << "pressure:" << endl;
for (int i = 0; i<100; i++) cout << this->pressure[i] << " " ;
vysledek.open("pressure" + to_string(step) + ".txt");
for (int i = 0; i<101; i++)
vysledek << 0.01*i << " " << pressure[i] << endl;
vysledek.close();
cout << " " << endl;
cout << "velocity:" << endl;
for (int i = 0; i<100; i++) cout << this->velocity[i] << " " ;
vysledek.open("velocity" + to_string(step) + ".txt");
for (int i = 0; i<101; i++)
vysledek << 0.01*i << " " << pressure[i] << endl;
vysledek.close();
cout << "energy:" << endl;
for (int i = 0; i<100; i++) cout << this->energy[i] << " " ;
vysledek.open("energy" + to_string(step) + ".txt");
for (int i = 0; i<101; i++)
vysledek << 0.01*i << " " << energy[i] << endl;
vysledek.close();
cout << " " << endl;
cout << "density:" << endl;
for (int i = 0; i<100; i++) cout << this->rho[i] << " " ;
vysledek.open("density" + to_string(step) + ".txt");
for (int i = 0; i<101; i++)
vysledek << 0.01*i << " " << rho[i] << endl;
vysledek.close();
getchar();
FileNameBaseNumberEnding( "rho-", step, 5, ".tnl", fileName );
if( ! rho.save( fileName ) )
return false;
FileNameBaseNumberEnding( "rhoVel-", step, 5, ".tnl", fileName );
if( ! rhoVel.save( fileName ) )
return false;
FileNameBaseNumberEnding( "energy-", step, 5, ".tnl", fileName );
if( ! energy.save( fileName ) )
return false;
return true;
}
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
void
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
getExplicitRHS( const RealType& time,
const RealType& tau,
const MeshType& mesh,
DofVectorType& _u,
DofVectorType& _fu,
MeshDependentDataType& meshDependentData )
{
/*
W[1] [0 ... count-1 ]
W[2] [count ... 2*count-1 ]
W[3] [2*count ... 3*count-1 ]
V this->velocity[]
p this->pressure[]
*/
typedef typename MeshType::Cell Cell;
int count = mesh.template getEntitiesCount< Cell >();
const RealType& size = 1;//mesh.template getSpaceStepsProducts< -1, 0 >();
for (long int i = 1; i < count - 1; i++)
_fu[i] = 1.0 / (2.0*tau) * (this->rho[i-1]+this->rho[i+1]-2.0*this->rho[i])-(1.0/(2.0 * size)) * (this->rho[i+1]
* this->velocity[i+1] - this->rho[i-1] * this->velocity[i-1]);
_fu[count-1] = _fu[count-2];
for (long int i = 1; i < count - 1; i++)
_fu[count + i] =
1.0 / (2.0 * tau) * (this->rhoVel[i+1] + this->rhoVel[i+1] -
2.0 * this->rhoVel[i])-(1.0 / (2.0 * size)) * ((this->rhoVel[i+1]
* this->velocity[i + 1] + this->pressure[i + 1]) - (this->rhoVel[i-1] * this->velocity[i - 1]
+ this->pressure[i - 1]));
_fu[2*count-1] = _fu[2*count-2];
for (long int i = 1; i < count - 1; i++)
_fu[i + 2 * count] = 1.0 / (2.0*tau) * (this->energy[i-1]
+ this->energy[i+1]-2.0*this->energy[i])
- (1.0/(2.0 * size)) * ((this->energy[i+1]
+ this->pressure[i + 1]) * this->velocity[i + 1] - (this->energy[i-1] + this->pressure[i -1]) * this->velocity[i - 1]);
_fu[3 * count-1] = _fu[3 * count-2];
for (long int i = 0; i <= count; i++) //pressure
this->pressure[i] = (this->gamma - 1 ) * ( this->energy[i] - 0.5 * this->rhoVel[i] * this->velocity[i]);
for (long int i = 0; i <= count ; i++) //velocity
this->velocity[i] = this->rhoVel[i]/this->rho[i];
/****
* If you use an explicit solver like tnlEulerSolver or tnlMersonSolver, you
* need to implement this method. Compute the right-hand side of
*
* d/dt u(x) = fu( x, u )
*
* You may use supporting mesh dependent data if you need.
this->bindDofs( mesh, _u );
tnlExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;
MeshFunctionType u( mesh, _u );
MeshFunctionType fu( mesh, _fu );
explicitUpdater.template update< typename Mesh::Cell >( time,
mesh,
this->differentialOperator,
this->boundaryCondition,
this->rightHandSide,
u,
fu );
tnlBoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter;
boundaryConditionsSetter.template apply< typename Mesh::Cell >(
this->boundaryCondition,
time + tau,
u ); */
}
template< typename Mesh,
typename BoundaryCondition,
typename RightHandSide,
typename DifferentialOperator >
template< typename Matrix >
void
eulerProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::
assemblyLinearSystem( const RealType& time,
const RealType& tau,
const MeshType& mesh,
DofVectorType& _u,
Matrix& matrix,
DofVectorType& b,
MeshDependentDataType& meshDependentData )
{
tnlLinearSystemAssembler< Mesh,
MeshFunctionType,
DifferentialOperator,
BoundaryCondition,
RightHandSide,
tnlBackwardTimeDiscretisation,
Matrix,
DofVectorType > systemAssembler;
tnlMeshFunction< Mesh > u( mesh, _u );
systemAssembler.template assembly< typename Mesh::Cell >( time,
tau,
mesh,
this->differentialOperator,
this->boundaryCondition,
this->rightHandSide,
u,
matrix,
b );
}
#endif /* eulerPROBLEM_IMPL_H_ */