Commit e0921e8d authored by Tomáš Oberhuber's avatar Tomáš Oberhuber
Browse files

Refactoring pressure getter.

parent 78cf27ea

examples/inviscid-flow/1d/EulerPressureGetter.h

+29 −184
Original line number Diff line number Diff line
@@ -3,24 +3,17 @@ 


#include <core/vectors/tnlVector.h>

#include <mesh/tnlGrid.h>

#include <functions/tnlDomain.h>



template< typename Mesh,

          typename Real = typename Mesh::RealType,

          typename Index = typename Mesh::IndexType >

class EulerPressureGetter

{

};



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

class EulerPressureGetter< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index >

: public tnlDomain< Mesh::geMeshDimensions(), MeshDomain >

{

   public:

      typedef tnlGrid< 1, MeshReal, Device, MeshIndex > MeshType;

      typedef typename MeshType::CoordinatesType CoordinatesType;

      

      typedef Mesh MeshType;

      typedef Real RealType;

      typedef Device DeviceType;

      typedef Index IndexType;
@@ -28,188 +21,40 @@ class EulerPressureGetter< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index 
      enum { Dimensions = MeshType::getMeshDimensions() };



      static tnlString getType();

      Real gamma;

      MeshFunctionType velocity;

      MeshFunctionType rhoVel;

      MeshFunctionType energy;

      

      void setGamma(const Real& gamma)

      {

          this->gamma = gamma;

      };



      void setVelocity(const MeshFunctionType& velocity)

      {

          this->velocity = velocity;

      };

      EulerPressureGetter( const MeshFunctionType& velocity,

                           const MeshFunctionType& rhoVel,

                           const MeshFunctionType& energy,

                           const RealType& gamma )

      : velocity( velocity ), rhoVel( rhoVel ), energy( energy ), gamma( gamma )

      {}



      void setRhoVel(const MeshFunctionType& rhoVel)

      {

          this->rhoVel = rhoVel;

      };



      void setEnergy(const MeshFunctionType& energy)

      {

          this->energy = energy;

      };



      template< typename MeshFunction, typename MeshEntity >

      __cuda_callable__

      Real operator()( const MeshFunction& u,

                       const MeshEntity& entity,

                       const RealType& time = 0.0 ) const;



      __cuda_callable__

      template< typename MeshEntity >

      Index getLinearSystemRowLength( const MeshType& mesh,

                                      const IndexType& index,

                                      const MeshEntity& entity ) const;



      template< typename MeshEntity, typename Vector, typename MatrixRow >

      __cuda_callable__

      void updateLinearSystem( const RealType& time,

                               const RealType& tau,

                               const MeshType& mesh,

                               const IndexType& index,

                               const MeshEntity& entity,

                               const MeshFunctionType& u,

                               Vector& b,

                               MatrixRow& matrixRow ) const;

};



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

class EulerPressureGetter< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index >

{

   public:

      typedef tnlGrid< 2, MeshReal, Device, MeshIndex > MeshType;

      typedef typename MeshType::CoordinatesType CoordinatesType;

      typedef Real RealType;

      typedef Device DeviceType;

      typedef Index IndexType;

      typedef tnlMeshFunction< MeshType > MeshFunctionType;

      enum { Dimensions = MeshType::getMeshDimensions() };



      static tnlString getType();

      Real gamma;

      MeshFunctionType velocity;

      MeshFunctionType rhoVel;

      MeshFunctionType energy;



      void setGamma(const Real& gamma)

      Real operator()( const MeshEntity& entity,

                       const RealType& time = 0.0 ) const

      {

          this->gamma = gamma;

      };



      void setVelocity(const MeshFunctionType& velocity)

      {

          this->velocity = velocity;

      };

         return this->operator[]( entity.getIndex() );

      }

      

      void setRhoVel(const MeshFunctionType& rhoVel)

      {

          this->rhoVel = rhoVel;

      };



      void setEnergy(const MeshFunctionType& energy)

      {

          this->energy = energy;

      };



      template< typename MeshFunction, typename MeshEntity >

      __cuda_callable__

      Real operator()( const MeshFunction& u,

                       const MeshEntity& entity,

                       const RealType& time = 0.0 ) const;



      __cuda_callable__

      template< typename MeshEntity >

      Index getLinearSystemRowLength( const MeshType& mesh,

                                      const IndexType& index,

                                      const MeshEntity& entity ) const;



      template< typename MeshEntity, typename Vector, typename MatrixRow >

      __cuda_callable__

      void updateLinearSystem( const RealType& time,

                               const RealType& tau,

                               const MeshType& mesh,

                               const IndexType& index,

                               const MeshEntity& entity,

                               const MeshFunctionType& u,

                               Vector& b,

                               MatrixRow& matrixRow ) const;

};



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

class EulerPressureGetter< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index >

      Real operator[]( const IndexType& idx ) const

      {

   public:

      typedef tnlGrid< 3, MeshReal, Device, MeshIndex > MeshType;

      typedef typename MeshType::CoordinatesType CoordinatesType;

      typedef Real RealType;

      typedef Device DeviceType;

      typedef Index IndexType;

      typedef tnlMeshFunction< MeshType > MeshFunctionType;

      enum { Dimensions = MeshType::getMeshDimensions() };

         return ( this->gamma - 1.0 ) * ( this->energy[ idx ] - 0.5 * this->rhoVel[ idx ] * this->velocity[ idx ]);

      }



      static tnlString getType();

      Real gamma;

      MeshFunctionType velocity;

      MeshFunctionType rhoVel;

      MeshFunctionType energy;

      

      void setGamma(const Real& gamma)

      {

          this->gamma = gamma;

      };

   protected:

      

      void setVelocity(const MeshFunctionType& velocity)

      {

          this->velocity = velocity;

      };



      void setRhoVel(const MeshFunctionType& rhoVel)

      {

          this->rhoVel = rhoVel;

      };

      Real gamma;

      

      void setEnergy(const MeshFunctionType& energy)

      {

          this->energy = energy;

      };

      const MeshFunctionType& velocity;

      

      template< typename MeshFunction, typename MeshEntity >

      __cuda_callable__

      Real operator()( const MeshFunction& u,

                       const MeshEntity& entity,

                       const RealType& time = 0.0 ) const;

      const MeshFunctionType& rhoVel;

      

      __cuda_callable__

      template< typename MeshEntity >

      Index getLinearSystemRowLength( const MeshType& mesh,

                                      const IndexType& index,

                                      const MeshEntity& entity ) const;

      const MeshFunctionType& energy;



      template< typename MeshEntity, typename Vector, typename MatrixRow >

      __cuda_callable__

      void updateLinearSystem( const RealType& time,

                               const RealType& tau,

                               const MeshType& mesh,

                               const IndexType& index,

                               const MeshEntity& entity,

                               const MeshFunctionType& u,

                               Vector& b,

                               MatrixRow& matrixRow ) const;

};





#include "EulerPressureGetter_impl.h"



#endif	/* EulerPressureGetter_H */

examples/inviscid-flow/1d/EulerPressureGetter_impl.h

deleted100644 → 0
+0 −336
Original line number Diff line number Diff line 
#ifndef EulerPressureGetter_IMPL_H

#define EulerPressureGetter_IMPL_H



/****

 * 1D problem

 */

template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

tnlString

EulerPressureGetter< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::

getType()

{

   return tnlString( "EulerPressureGetter< " ) +

          MeshType::getType() + ", " +

          ::getType< Real >() + ", " +

          ::getType< Index >() + " >";

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

template< typename MeshFunction, typename MeshEntity >

__cuda_callable__

Real

EulerPressureGetter< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::

operator()( const MeshFunction& u,

            const MeshEntity& entity,

            const Real& time ) const

{

    static_assert( MeshEntity::entityDimensions == 1, "Wrong mesh entity dimensions." ); 

    static_assert( MeshFunction::getEntitiesDimensions() == 1, "Wrong preimage function" ); 

    const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities();   

   //pressure

   const IndexType& center = entity.getIndex(); 

   return ( this->gamma - 1 ) * ( energy[ center ] - 0.5 * rhoVel[ center ] * velocity[ center ]);

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

template< typename MeshEntity >

__cuda_callable__

Index

EulerPressureGetter< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::

getLinearSystemRowLength( const MeshType& mesh,

                          const IndexType& index,

                          const MeshEntity& entity ) const

{

   /****

    * Return a number of non-zero elements in a line (associated with given grid element) of

    * the linear system.

    * The following example is the Laplace operator approximated 

    * by the Finite difference method.

    */



   return 2*Dimensions + 1;

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

   template< typename MeshEntity, typename Vector, typename MatrixRow >

__cuda_callable__

void

EulerPressureGetter< tnlGrid< 1, MeshReal, Device, MeshIndex >, Real, Index >::

updateLinearSystem( const RealType& time,

                    const RealType& tau,

                    const MeshType& mesh,

                    const IndexType& index,

                    const MeshEntity& entity,

                    const MeshFunctionType& u,

                    Vector& b,

                    MatrixRow& matrixRow ) const

{

   /****

    * Setup the non-zero elements of the linear system here.

    * The following example is the Laplace operator appriximated 

    * by the Finite difference method.

    */



    const typename MeshEntity::template NeighbourEntities< 1 >& neighbourEntities = entity.getNeighbourEntities(); 

   const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2 >(); 

   const IndexType& center = entity.getIndex(); 

   const IndexType& east = neighbourEntities.template getEntityIndex< 1 >(); 

   const IndexType& west = neighbourEntities.template getEntityIndex< -1 >(); 

   matrixRow.setElement( 0, west,   - lambdaX );

   matrixRow.setElement( 1, center, 2.0 * lambdaX );

   matrixRow.setElement( 2, east,   - lambdaX );

}



/****

 * 2D problem

 */

template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

tnlString

EulerPressureGetter< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::

getType()

{

   return tnlString( "EulerPressureGetter< " ) +

          MeshType::getType() + ", " +

          ::getType< Real >() + ", " +

          ::getType< Index >() + " >";

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

template< typename MeshFunction, typename MeshEntity >

__cuda_callable__

Real

EulerPressureGetter< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::

operator()( const MeshFunction& u,

            const MeshEntity& entity,

            const Real& time ) const

{

   /****

    * Implement your explicit form of the differential operator here.

    * The following example is the Laplace operator approximated 

    * by the Finite difference method.

    */

    static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); 

    static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); 

    const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); 



   const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); 

   const RealType& hySquareInverse = entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); 

   const IndexType& center = entity.getIndex(); 

   const IndexType& east  = neighbourEntities.template getEntityIndex<  1,  0 >(); 

   const IndexType& west  = neighbourEntities.template getEntityIndex< -1,  0 >(); 

   const IndexType& north = neighbourEntities.template getEntityIndex<  0,  1 >(); 

   const IndexType& south = neighbourEntities.template getEntityIndex<  0, -1 >(); 

   return ( u[ west ] - 2.0 * u[ center ] + u[ east ]  ) * hxSquareInverse +

          ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse;

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

template< typename MeshEntity >

__cuda_callable__

Index

EulerPressureGetter< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::

getLinearSystemRowLength( const MeshType& mesh,

                          const IndexType& index,

                          const MeshEntity& entity ) const

{

   /****

    * Return a number of non-zero elements in a line (associated with given grid element) of

    * the linear system.

    * The following example is the Laplace operator approximated 

    * by the Finite difference method.

    */



   return 2*Dimensions + 1;

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

   template< typename MeshEntity, typename Vector, typename MatrixRow >

__cuda_callable__

void

EulerPressureGetter< tnlGrid< 2, MeshReal, Device, MeshIndex >, Real, Index >::

updateLinearSystem( const RealType& time,

                    const RealType& tau,

                    const MeshType& mesh,

                    const IndexType& index,

                    const MeshEntity& entity,

                    const MeshFunctionType& u,

                    Vector& b,

                    MatrixRow& matrixRow ) const

{

   /****

    * Setup the non-zero elements of the linear system here.

    * The following example is the Laplace operator appriximated 

    * by the Finite difference method.

    */



    const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); 

   const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2, 0 >(); 

   const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts< 0, -2 >(); 

   const IndexType& center = entity.getIndex(); 

   const IndexType& east  = neighbourEntities.template getEntityIndex<  1,  0 >(); 

   const IndexType& west  = neighbourEntities.template getEntityIndex< -1,  0 >(); 

   const IndexType& north = neighbourEntities.template getEntityIndex<  0,  1 >(); 

   const IndexType& south = neighbourEntities.template getEntityIndex<  0, -1 >(); 

   matrixRow.setElement( 0, south,  -lambdaY );

   matrixRow.setElement( 1, west,   -lambdaX );

   matrixRow.setElement( 2, center, 2.0 * ( lambdaX + lambdaY ) );

   matrixRow.setElement( 3, east,   -lambdaX );

   matrixRow.setElement( 4, north,  -lambdaY );

}



/****

 * 3D problem

 */

template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

tnlString

EulerPressureGetter< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::

getType()

{

   return tnlString( "EulerPressureGetter< " ) +

          MeshType::getType() + ", " +

          ::getType< Real >() + ", " +

          ::getType< Index >() + " >";

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

template< typename MeshFunction, typename MeshEntity >

__cuda_callable__

Real

EulerPressureGetter< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::

operator()( const MeshFunction& u,

            const MeshEntity& entity,

            const Real& time ) const

{

   /****

    * Implement your explicit form of the differential operator here.

    * The following example is the Laplace operator approximated 

    * by the Finite difference method.

    */

    static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); 

    static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); 

    const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); 



   const RealType& hxSquareInverse = entity.getMesh().template getSpaceStepsProducts< -2,  0,  0 >(); 

   const RealType& hySquareInverse = entity.getMesh().template getSpaceStepsProducts<  0, -2,  0 >(); 

   const RealType& hzSquareInverse = entity.getMesh().template getSpaceStepsProducts<  0,  0, -2 >(); 

   const IndexType& center = entity.getIndex(); 

   const IndexType& east  = neighbourEntities.template getEntityIndex<  1,  0,  0 >(); 

   const IndexType& west  = neighbourEntities.template getEntityIndex< -1,  0,  0 >(); 

   const IndexType& north = neighbourEntities.template getEntityIndex<  0,  1,  0 >(); 

   const IndexType& south = neighbourEntities.template getEntityIndex<  0, -1,  0 >(); 

   const IndexType& up    = neighbourEntities.template getEntityIndex<  0,  0,  1 >(); 

   const IndexType& down  = neighbourEntities.template getEntityIndex<  0,  0, -1 >(); 

   return ( u[ west ] - 2.0 * u[ center ] + u[ east ]  ) * hxSquareInverse +

          ( u[ south ] - 2.0 * u[ center ] + u[ north ] ) * hySquareInverse +

          ( u[ up ] - 2.0 * u[ center ] + u[ down ] ) * hzSquareInverse;

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

template< typename MeshEntity >

__cuda_callable__

Index

EulerPressureGetter< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::

getLinearSystemRowLength( const MeshType& mesh,

                          const IndexType& index,

                          const MeshEntity& entity ) const

{

   /****

    * Return a number of non-zero elements in a line (associated with given grid element) of

    * the linear system.

    * The following example is the Laplace operator approximated 

    * by the Finite difference method.

    */



   return 2*Dimensions + 1;

}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

   template< typename MeshEntity, typename Vector, typename MatrixRow >

__cuda_callable__

void

EulerPressureGetter< tnlGrid< 3, MeshReal, Device, MeshIndex >, Real, Index >::

updateLinearSystem( const RealType& time,

                    const RealType& tau,

                    const MeshType& mesh,

                    const IndexType& index,

                    const MeshEntity& entity,

                    const MeshFunctionType& u,

                    Vector& b,

                    MatrixRow& matrixRow ) const

{

   /****

    * Setup the non-zero elements of the linear system here.

    * The following example is the Laplace operator appriximated 

    * by the Finite difference method.

    */



    const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); 

   const RealType& lambdaX = tau * entity.getMesh().template getSpaceStepsProducts< -2,  0,  0 >(); 

   const RealType& lambdaY = tau * entity.getMesh().template getSpaceStepsProducts<  0, -2,  0 >(); 

   const RealType& lambdaZ = tau * entity.getMesh().template getSpaceStepsProducts<  0,  0, -2 >(); 

   const IndexType& center = entity.getIndex(); 

   const IndexType& east  = neighbourEntities.template getEntityIndex<  1,  0,  0 >(); 

   const IndexType& west  = neighbourEntities.template getEntityIndex< -1,  0,  0 >(); 

   const IndexType& north = neighbourEntities.template getEntityIndex<  0,  1,  0 >(); 

   const IndexType& south = neighbourEntities.template getEntityIndex<  0, -1,  0 >(); 

   const IndexType& up    = neighbourEntities.template getEntityIndex<  0,  0,  1 >(); 

   const IndexType& down  = neighbourEntities.template getEntityIndex<  0,  0, -1 >(); 

   matrixRow.setElement( 0, down,   -lambdaZ );

   matrixRow.setElement( 1, south,  -lambdaY );

   matrixRow.setElement( 2, west,   -lambdaX );

   matrixRow.setElement( 3, center, 2.0 * ( lambdaX + lambdaY + lambdaZ ) );

   matrixRow.setElement( 4, east,   -lambdaX );

   matrixRow.setElement( 5, north,  -lambdaY );

   matrixRow.setElement( 6, up,     -lambdaZ );

}



#endif	/* EulerPressureGetterIMPL_H */

examples/inviscid-flow/1d/eulerProblem.h

+9 −15
Original line number Diff line number Diff line
@@ -33,24 +33,9 @@ class eulerProblem: 
      typedef typename DifferentialOperator::Velocity Velocity;

      typedef typename DifferentialOperator::Pressure Pressure;

      

    //definition

	tnlVector< RealType, DeviceType, IndexType > _uRho;

	tnlVector< RealType, DeviceType, IndexType > _uRhoVelocity;

	tnlVector< RealType, DeviceType, IndexType > _uEnergy;



	tnlVector< RealType, DeviceType, IndexType > _fuRho;

	tnlVector< RealType, DeviceType, IndexType > _fuRhoVelocity;

	tnlVector< RealType, DeviceType, IndexType > _fuEnergy;





      static tnlString getTypeStatic();

      tnlVector< RealType, DeviceType, IndexType > data;

      tnlVector< RealType, DeviceType, IndexType > rho;

      tnlVector< RealType, DeviceType, IndexType > rhoVel;

      tnlVector< RealType, DeviceType, IndexType > energy;

      tnlVector< RealType, DeviceType, IndexType > pressure;

      tnlVector< RealType, DeviceType, IndexType > velocity;

      double gamma;



      tnlString getPrologHeader() const;
@@ -106,6 +91,15 @@ class eulerProblem: 
      DifferentialOperator differentialOperator;

      BoundaryCondition boundaryCondition;

      RightHandSide rightHandSide;

      

      MeshFunctionType uRho, uRhoVelocity, uEnergy;

      MeshFunctionType fuRho, fuRhoVelocity, fuEnergy;

      

      MeshFunctionType rho, rhoVel, energy, pressure, velocity;

      

      RealType gamma;





};



#include "eulerProblem_impl.h"

examples/inviscid-flow/1d/eulerProblem_impl.h

+16 −16
Original line number Diff line number Diff line
@@ -246,30 +246,30 @@ getExplicitRHS( const RealType& time, 
                MeshDependentDataType& meshDependentData )

{

    typedef typename MeshType::Cell Cell;

    int count = mesh.template getEntitiesCount< Cell >()/3;

    int count = mesh.template getEntitiesCount< Cell >(); ///3;

	//bind _u

    this->_uRho.bind(_u,0,count);

    this->_uRhoVelocity.bind(_u,count,count);

    this->_uEnergy.bind(_u,2 * count,count);

    this->uRho.bind(_u,0,count);

    this->uRhoVelocity.bind(_u,count,count);

    this->uEnergy.bind(_u,2 * count,count);

		

	//bind _fu

    this->_fuRho.bind(_u,0,count);

    this->_fuRhoVelocity.bind(_u,count,count);

    this->_fuEnergy.bind(_u,2 * count,count);

    this->fuRho.bind(_u,0,count);

    this->fuRhoVelocity.bind(_u,count,count);

    this->fuEnergy.bind(_u,2 * count,count);

        //bind MeshFunctions

   MeshFunctionType velocity( mesh, this->velocity );

   MeshFunctionType pressure( mesh, this->pressure );

   MeshFunctionType uRho( mesh, _uRho ); 

   MeshFunctionType fuRho( mesh, _fuRho );

   MeshFunctionType uRhoVelocity( mesh, _uRhoVelocity ); 

   MeshFunctionType fuRhoVelocity( mesh, _fuRhoVelocity );

   MeshFunctionType uEnergy( mesh, _uEnergy ); 

   MeshFunctionType fuEnergy( mesh, _fuEnergy );

   this->velocity.setMesh( mesh );

   // TODO: osetrit rychlost podobne jako tlak

   this->pressure.setMesh( mesh );

   EulerPressureGetter< Mesh, RealType, IndexType > pressureGetter( velocity, rhoVel, energy, gamma );

   this->pressure = pressureGetter;

   //generating Differential operator object

   Continuity lF1DContinuity;

   Momentum lF1DMomentum;

   Energy lF1DEnergy;



   

   

   

   //rho

   this->bindDofs( mesh, _u );

   lF1DContinuity.setTau(tau);