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

Fixing the inviscid flow solver.

parent 13bde21f
Loading
Loading
Loading
Loading
+1 −1
Original line number Diff line number Diff line
@@ -6,7 +6,7 @@ set( tnl_inviscid_flow_SOURCES
     euler.cu )
               
IF( BUILD_CUDA )
   CUDA_ADD_EXECUTABLE(tnl-euler-2d${debugExt} euler-cuda.cu)
   CUDA_ADD_EXECUTABLE(tnl-euler-2d${debugExt} euler.cu)
   target_link_libraries (tnl-euler-2d${debugExt} tnl${debugExt}-${tnlVersion}  ${CUSPARSE_LIBRARY} )
ELSE(  BUILD_CUDA )               
   ADD_EXECUTABLE(tnl-euler-2d${debugExt} euler.cpp)     
+1 −2
Original line number Diff line number Diff line
@@ -140,9 +140,8 @@ class CompressibleConservativeVariables
      
      MeshFunctionPointer density;
      MomentumFieldPointer momentum;
      //MeshFunctionPointer pressure;
      MeshFunctionPointer energy;
      
};

} // namespace TNL
} // namespace TN
 No newline at end of file
+9 −9
Original line number Diff line number Diff line
@@ -109,8 +109,8 @@ class LaxFridrichsContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Re
         const IndexType& west = neighbourEntities.template getEntityIndex< -1 >();
         const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
         const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
         return (0.5 * this->tau) * ( u[ west ] - 2.0 * u[ center ]  + u[ east ] ) 
               - 0.5 * hxInverse * ( u[ west ] * velocity_x_west - u[ east ] * velocity_x_east );
         return 1.0 / 2.0 * this->tau * ( u[ west ] - 2.0 * u[ center ]  + u[ east ] ) 
               - 0.5 * ( u[ west ] * velocity_x_west - u[ east ] * velocity_x_east ) * hxInverse;
      }

      /*template< typename MeshEntity >
@@ -175,9 +175,9 @@ class LaxFridrichsContinuity< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Re
         const RealType& velocity_y_north = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ];
         const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ];
         
         return 0.5 * ( this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] ) 
                       - ( u[ west ] * velocity_x_west - u[ east ] * velocity_x_east ) * hxInverse
                       - ( u[ north ] * velocity_y_north - u[ south ] * velocity_y_south ) * hyInverse );
         return 1.0 / 4.0 * this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] ) 
                       - 0.5 * ( ( u[ west ] * velocity_x_west - u[ east ] * velocity_x_east ) * hxInverse
                               + ( u[ north ] * velocity_y_north - u[ south ] * velocity_y_south ) * hyInverse );
      }

      /*template< typename MeshEntity >
@@ -248,10 +248,10 @@ class LaxFridrichsContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Re
         const RealType& velocity_z_up    = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ];
         const RealType& velocity_z_down  = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ];
         
         return 0.5 * ( this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] + u[ up ] + u[ down ]- 6.0 * u[ center ] ) 
                        - ( u[ west ] * velocity_x_west - u[ east ] * velocity_x_east ) * hxInverse
                        - ( u[ north ] * velocity_y_north - u[ south ] * velocity_y_south ) * hyInverse
                        - ( u[ up ] * velocity_z_up - u[ down ] * velocity_z_down ) * hzInverse );
         return 1.0 / 6.0 * this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] + u[ up ] + u[ down ]- 6.0 * u[ center ] ) 
                - 0.5 * ( ( u[ west ] * velocity_x_west - u[ east ] * velocity_x_east ) * hxInverse
                        + ( u[ north ] * velocity_y_north - u[ south ] * velocity_y_south ) * hyInverse
                        + ( u[ up ] * velocity_z_up - u[ down ] * velocity_z_down ) * hzInverse );
         
      }

+13 −13
Original line number Diff line number Diff line
@@ -113,9 +113,9 @@ class LaxFridrichsEnergy< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real,
         const RealType& pressure_east = this->pressure.template getData< DeviceType >()[ east ];
         const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
         const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
         return 0.5 * ( this->tau * ( e[ west ] - 2.0 * e[ center ]  + e[ east ] ) 
                       - ( ( e[ west ] + pressure_west ) * velocity_x_west  
                         - ( e[ east ] + pressure_east ) * velocity_x_east ) * hxInverse );
         return 1.0 / 2.0 * this->tau * ( e[ west ] - 2.0 * e[ center ]  + e[ east ] ) 
                - 0.5 * ( ( e[ west ] + pressure_west ) * velocity_x_west  
                         - ( e[ east ] + pressure_east ) * velocity_x_east ) * hxInverse;
         
      }

@@ -186,10 +186,10 @@ class LaxFridrichsEnergy< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real,
         const RealType& velocity_y_north = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ];
         const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ];         
         
         return 0.5 * ( this->tau * ( e[ west ] + e[ east ] + e[ south ] + e[ north ] - 4.0 * e[ center ] ) 
                       - ( ( ( e[ west ] + pressure_west ) * velocity_x_west )
         return 1.0 / 4.0 * this->tau * ( e[ west ] + e[ east ] + e[ south ] + e[ north ] - 4.0 * e[ center ] ) 
                - 0.5 * ( ( ( ( e[ west ] + pressure_west ) * velocity_x_west )
                          -( ( e[ east ] + pressure_east ) * velocity_x_east ) ) * hxInverse
                       - ( ( ( e[ north ] + pressure_north ) * velocity_y_north )
                        + ( ( ( e[ north ] + pressure_north ) * velocity_y_north )
                          -( ( e[ south ] + pressure_south ) * velocity_y_south ) ) * hyInverse );
      }

@@ -268,12 +268,12 @@ class LaxFridrichsEnergy< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real,
         const RealType& velocity_z_up    = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ];
         const RealType& velocity_z_down  = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ];         
         
         return 0.5 * ( this->tau * ( e[ west ] + e[ east ] + e[ south ] + e[ north ] + e[ up ] + e[ down ] - 6.0 * e[ center ] ) 
                       - ( ( ( e[ west ] + pressure_west ) * velocity_x_west )
         return 1.0 / 6.0 * this->tau * ( e[ west ] + e[ east ] + e[ south ] + e[ north ] + e[ up ] + e[ down ] - 6.0 * e[ center ] ) 
                - 0.5 * ( ( ( ( e[ west ] + pressure_west ) * velocity_x_west )
                           -( ( e[ east ] + pressure_east ) * velocity_x_east ) ) * hxInverse
                       - ( ( ( e[ north ] + pressure_north ) * velocity_y_north )
                        + ( ( ( e[ north ] + pressure_north ) * velocity_y_north )
                           -( ( e[ south ] + pressure_south ) * velocity_y_south ) ) * hyInverse
                       - ( ( ( e[ up ] + pressure_up ) * velocity_z_up )
                        + ( ( ( e[ up ] + pressure_up ) * velocity_z_up )
                           -( ( e[ down ] + pressure_down ) * velocity_z_down ) ) * hzInverse );
      }

+18 −18
Original line number Diff line number Diff line
@@ -58,7 +58,7 @@ class LaxFridrichsMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Rea

      template< typename MeshFunction, typename MeshEntity >
      __cuda_callable__
      Real operator()( const MeshFunction& u,
      Real operator()( const MeshFunction& rho_u,
                       const MeshEntity& entity,
                       const RealType& time = 0.0 ) const
      {
@@ -75,9 +75,9 @@ class LaxFridrichsMomentumX< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Rea
         const RealType& velocity_x_east = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];
         const RealType& velocity_x_west = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ west ];
         
         return 0.5 * ( this->tau * ( u[ west ]  + u[ east ]  - 2.0 * u[ center ] ) 
                        - ( ( u[ west ] * velocity_x_west + pressure_west ) 
                          - ( u[ east ] * velocity_x_east + pressure_east ) ) * hxInverse );
         return 1.0 / 2.0 * this->tau * ( rho_u[ west ]  + rho_u[ east ]  - 2.0 * rho_u[ center ] ) 
                - 0.5 * ( ( rho_u[ west ] * velocity_x_west + pressure_west ) 
                         -( rho_u[ east ] * velocity_x_east + pressure_east ) ) * hxInverse;
      }

      /*template< typename MeshEntity >
@@ -130,7 +130,7 @@ class LaxFridrichsMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Rea

      template< typename MeshFunction, typename MeshEntity >
      __cuda_callable__
      Real operator()( const MeshFunction& u,
      Real operator()( const MeshFunction& rho_u,
                       const MeshEntity& entity,
                       const RealType& time = 0.0 ) const
      {
@@ -154,11 +154,11 @@ class LaxFridrichsMomentumX< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Rea
         const RealType& velocity_y_north = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ north ];
         const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ];         
         
         return 0.5 * ( this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4.0 * u[ center ] ) 
                       - ( ( u[ west ] * velocity_x_west + pressure_west )
                         - ( u[ east ] * velocity_x_east + pressure_east ) )* hxInverse
                       - ( ( u[ north ] * velocity_y_north )
                         - ( u[ south ] * velocity_y_south ) )* hyInverse );
         return 1.0 / 4.0 * this->tau * ( rho_u[ west ] + rho_u[ east ] + rho_u[ south ] + rho_u[ north ] - 4.0 * rho_u[ center ] ) 
                - 0.5 * ( ( ( rho_u[ west ] * velocity_x_west + pressure_west )
                          - ( rho_u[ east ] * velocity_x_east + pressure_east ) ) * hxInverse
                        + ( ( rho_u[ north ] * velocity_y_north )
                          - ( rho_u[ south ] * velocity_y_south ) ) * hyInverse );
      }

      /*template< typename MeshEntity >
@@ -211,7 +211,7 @@ class LaxFridrichsMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real

      template< typename MeshFunction, typename MeshEntity >
      __cuda_callable__
      Real operator()( const MeshFunction& u,
      Real operator()( const MeshFunction& rho_u,
                       const MeshEntity& entity,
                       const RealType& time = 0.0 ) const
      {
@@ -243,13 +243,13 @@ class LaxFridrichsMomentumX< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real
         const RealType& velocity_y_south = this->velocity.template getData< DeviceType >()[ 1 ].template getData< DeviceType >()[ south ];
         const RealType& velocity_z_up    = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ up ];
         const RealType& velocity_z_down  = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ];
         return 0.5 * ( this->tau * ( u[ west ] + u[ east ] + u[ south ] + u[ north ] + u[ up ] + u[ down ] - 6.0 * u[ center ] ) 
                       - ( ( u[ west ] * velocity_x_west + pressure_west )
                         - ( u[ east ] * velocity_x_east + pressure_east ) )* hxInverse
                       - ( ( u[ north ] * velocity_y_north )
                         - ( u[ south ] * velocity_y_south ) )* hyInverse
                       - ( ( u[ up ] * velocity_z_up )
                         - ( u[ down ] * velocity_z_down ) )* hzInverse );
         return 1.0 / 6.0 * this->tau * ( rho_u[ west ] + rho_u[ east ] + rho_u[ south ] + rho_u[ north ] + rho_u[ up ] + rho_u[ down ] - 6.0 * rho_u[ center ] ) 
                - 0.5 * ( ( ( rho_u[ west ] * velocity_x_west + pressure_west )
                          - ( rho_u[ east ] * velocity_x_east + pressure_east ) )* hxInverse
                        + ( ( rho_u[ north ] * velocity_y_north )
                          - ( rho_u[ south ] * velocity_y_south ) )* hyInverse
                        + ( ( rho_u[ up ] * velocity_z_up )
                          - ( rho_u[ down ] * velocity_z_down ) )* hzInverse );
      }

      /*template< typename MeshEntity >
Loading