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

Adding inviscid flow and advection.

parent 6b078884

examples/CMakeLists.txt

+2 −0
Original line number Diff line number Diff line 
add_subdirectory( heat-equation )

add_subdirectory( navier-stokes )

add_subdirectory( advection )

add_subdirectory( inviscid-flow )

#add_subdirectory( mean-curvature-flow )

examples/advection/CMakeLists.txt

0 → 100644
+20 −0
Original line number Diff line number Diff line 
set( tnl_heat_equation_SOURCES     

     advection.cpp

     advection.cu )

               

IF( BUILD_CUDA )

   CUDA_ADD_EXECUTABLE( tnl-advection{debugExt} advection.cu)   

   target_link_libraries( tnl-advection${debugExt} tnl${debugExt}-${tnlVersion}  ${CUSPARSE_LIBRARY} )

ELSE(  BUILD_CUDA )               

   ADD_EXECUTABLE( tnl-advection${debugExt} advection.cpp)     

   target_link_libraries( tnl-advection${debugExt} tnl${debugExt}-${tnlVersion} )

ENDIF( BUILD_CUDA )





INSTALL( TARGETS tnl-advection${debugExt}

         RUNTIME DESTINATION bin

         PERMISSIONS OWNER_READ OWNER_WRITE OWNER_EXECUTE GROUP_READ GROUP_EXECUTE WORLD_READ WORLD_EXECUTE )

        

INSTALL( FILES tnl-run-advection

               ${tnl_heat_equation_SOURCES}

         DESTINATION share/tnl-${tnlVersion}/examples/advection )

examples/advection/LaxFridrichs.h

0 → 100644
+218 −0
Original line number Diff line number Diff line 
#ifndef LaxFridrichs_H

#define LaxFridrichs_H



#include <core/vectors/tnlVector.h>

#include <mesh/tnlGrid.h>



template< typename Mesh,

          typename Real = typename Mesh::RealType,

          typename Index = typename Mesh::IndexType >

class LaxFridrichs

{

};



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

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

{

   public:

      typedef tnlGrid< 1, 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() };

      Real tau;

      double artificalViscosity;

      double advectionSpeedX;

      double advectionSpeedY;



      void setAdvectionSpeedY(const double advectionSpeed)

      {

	   this->advectionSpeedY = advectionSpeed;

      }





      void setAdvectionSpeedX(const double advectionSpeed)

      {

	   this->advectionSpeedX = advectionSpeed;

      }



      void setViscosity(const double artificalViscosity)

      {

	   this->artificalViscosity = artificalViscosity;

      }

      

      void setTau(const Real& tau)

      {

          this->tau = tau;

      };



      static tnlString getType();



      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 LaxFridrichs< 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() };

      Real tau;

      double artificalViscosity;

      double advectionSpeedX;

      double advectionSpeedY;



      void setAdvectionSpeedY(const double advectionSpeed)

      {

	   this->advectionSpeedY = advectionSpeed;

      }





      void setAdvectionSpeedX(const double advectionSpeed)

      {

	   this->advectionSpeedX = advectionSpeed;

      }



      void setViscosity(const double artificalViscosity)

      {

	   this->artificalViscosity = artificalViscosity;

      }

      

      void setTau(const Real& tau)

      {

          this->tau = tau;

      };



      static tnlString getType();



      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 LaxFridrichs< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index >

{

   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() };

      Real tau;

      double artificalViscosity;

      double advectionSpeedX;

      double advectionSpeedY;



      void setAdvectionSpeedY(const double advectionSpeed)

      {

	   this->advectionSpeedY = advectionSpeed;

      }





      void setAdvectionSpeedX(const double advectionSpeed)

      {

	   this->advectionSpeedX = advectionSpeed;

      }



      void setViscosity(const double artificalViscosity)

      {

	   this->artificalViscosity = artificalViscosity;

      }

      

      void setTau(const Real& tau)

      {

          this->tau = tau;

      };



      static tnlString getType();



      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;

};





#include "LaxFridrichs_impl.h"



#endif	/* LaxFridrichs_H */

examples/advection/LaxFridrichs_impl.h

0 → 100644
+369 −0
Original line number Diff line number Diff line 
#ifndef LaxFridrichs_IMPL_H

#define LaxFridrichs_IMPL_H



/****

 * 1D problem

 */

template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

tnlString

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

getType()

{

   return tnlString( "LaxFridrichs< " ) +

          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

LaxFridrichs< tnlGrid< 1, 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 == 1, "Wrong mesh entity dimensions." ); 

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

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



   const RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1 >(); 

   const IndexType& center = entity.getIndex(); 

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

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

   double a;

   a = this->advectionSpeedX;

   return   (0.5 / this->tau ) * this->artificalViscosity *

	    ( u[ west ]

            - 2.0 * u[ center ]

            + u[ east ] )

            - (a * ( u[ east ]

            - u[west] ) )

            * hxInverse * 0.5;



}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

template< typename MeshEntity >

__cuda_callable__

Index

LaxFridrichs< 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

LaxFridrichs< 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

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

getType()

{

   return tnlString( "LaxFridrichs< " ) +

          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

LaxFridrichs< 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& hxInverse = entity.getMesh().template getSpaceStepsProducts< -1, 0 >(); 

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

   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 >(); 

   double a;

   double b;

   a = this->advectionSpeedX;

   b = this->advectionSpeedY;

   return ( 0.25 / this->tau ) * this->artificalViscosity * (

              u[ west ]

            + u[ east ]

            + u[ south ]

            + u[ north ]

            - 4 * u[ center ] )

            - (a * ( u[ east ]

            - u[west] ) )

            * hxInverse * 0.5

            - (b * ( u[ south ]

            - u[ south ] ) )

            * hyInverse * 0.5;





}



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

template< typename MeshEntity >

__cuda_callable__

Index

LaxFridrichs< 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

LaxFridrichs< 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

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

getType()

{

   return tnlString( "LaxFridrichs< " ) +

          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

LaxFridrichs< 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

LaxFridrichs< 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

LaxFridrichs< 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	/* LaxFridrichsIMPL_H */

examples/advection/advection.cpp

0 → 100644
+1 −0
Original line number Diff line number Diff line 
#include "advection.h"