Commit 6093dd3c authored by Jan Schäfer's avatar Jan Schäfer
Browse files

Added 3d computations for stegger-warming, cleared unused velocity binding in... 

Added 3d computations for stegger-warming, cleared unused velocity binding in stegger-warming algoryrhm, specified energy computations for each dimension
parent dfd39fbf

examples/CMakeLists.txt

+1 −0
Original line number Diff line number Diff line
@@ -14,3 +14,4 @@ add_subdirectory( inviscid-flow ) 
add_subdirectory( inviscid-flow-sw )

#add_subdirectory( mean-curvature-flow )

add_subdirectory( flow )

add_subdirectory( flow-sw )

examples/inviscid-flow-sw/Upwind.h

+8 −0
Original line number Diff line number Diff line
@@ -95,6 +95,14 @@ class Upwind 
         this->energyOperatorPointer->setPressure( pressure );

      }



      void setDensity( const MeshFunctionPointer& density )

      {

         this->momentumXOperatorPointer->setDensity( density );

         this->momentumYOperatorPointer->setDensity( density );

         this->momentumZOperatorPointer->setDensity( density );

         this->energyOperatorPointer->setDensity( density );

      }

      

      void setVelocity( const VectorFieldPointer& velocity )

      {

         this->continuityOperatorPointer->setVelocity( velocity );

examples/inviscid-flow-sw/UpwindContinuity.h

+15 −15
Original line number Diff line number Diff line
@@ -344,22 +344,22 @@ class UpwindContinuity< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, In 
         const RealType& velocity_z_down   = this->velocity.template getData< DeviceType >()[ 2 ].template getData< DeviceType >()[ down ];

         

         return -hxInverse * (

                                  u[ center ] / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )

                                - u[ west ]   / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_x_west   + std::sqrt( this->gamma * pressure_west   / u[ west ]   ) )

                                - u[ center ] / ( 2 * this->gamma ) * ( velocity_x_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )

                                + u[ east ]   / ( 2 * this->gamma ) * ( velocity_x_east   - std::sqrt( this->gamma * pressure_east   / u[ east   ] ) )

                                   this->positiveDensityFlux( u[ center ], velocity_x_center, pressure_center )

                                -  this->positiveDensityFlux( u[ west   ], velocity_x_west  , pressure_west   )

                                -  this->negativeDensityFlux( u[ center ], velocity_x_center, pressure_center )

                                +  this->negativeDensityFlux( u[ east   ], velocity_x_east  , pressure_east   )

                             )

                -hyInverse * (

                                  u[ center ] / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_y_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )

                                - u[ south ]  / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_y_south  + std::sqrt( this->gamma * pressure_south  / u[ south ]  ) )

                                - u[ center ] / ( 2 * this->gamma ) * ( velocity_y_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) )

                                + u[ north ]  / ( 2 * this->gamma ) * ( velocity_y_north  - std::sqrt( this->gamma * pressure_north  / u[ north ]  ) )

                                   this->positiveDensityFlux( u[ center ], velocity_y_center, pressure_center )

                                -  this->positiveDensityFlux( u[ south  ], velocity_y_south , pressure_south  )

                                -  this->negativeDensityFlux( u[ center ], velocity_y_center, pressure_center )

                                +  this->negativeDensityFlux( u[ north  ], velocity_y_north , pressure_north  )

                             )

                -hzInverse * (

                                  u[ center ] / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_z_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )

                                - u[ down ]   / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_z_down   + std::sqrt( this->gamma * pressure_down   / u[ down ]   ) )

                                - u[ center ] / ( 2 * this->gamma ) * ( velocity_z_center - std::sqrt( this->gamma * pressure_center / u[ center ] ) ) 

                                + u[ up ]     / ( 2 * this->gamma ) * ( velocity_z_up     - std::sqrt( this->gamma * pressure_up     / u[ up ]     ) ) 

                                   this->positiveDensityFlux( u[ center ], velocity_z_center, pressure_center )

                                -  this->positiveDensityFlux( u[ down   ], velocity_z_down  , pressure_down   )

                                -  this->negativeDensityFlux( u[ center ], velocity_z_center, pressure_center )

                                +  this->negativeDensityFlux( u[ up     ], velocity_z_up    , pressure_up     )

                             );

         

      }

examples/inviscid-flow-sw/UpwindEnergy.h

+197 −186

File changed.

Preview size limit exceeded, changes collapsed.

examples/inviscid-flow-sw/UpwindMomentumBase.h

+11 −4
Original line number Diff line number Diff line
@@ -47,6 +47,11 @@ class UpwindMomentumBase 
          this->velocity = velocity;

      };



      void setDensity( const MeshFunctionPointer& density )

      {

          this->density = density;

      };

      

      void setPressure( const MeshFunctionPointer& pressure )

      {

          this->pressure = pressure;
@@ -85,7 +90,7 @@ class UpwindMomentumBase 
         const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );

         const RealType& machNumber = velocity_main / speedOfSound;

         if ( machNumber <= -1.0 )

            return 0;

            return 0.0;

        else if ( machNumber <= 0.0 )

            return density * speedOfSound / ( 2 * this->gamma ) * ( machNumber + 1.0 ) * velocity_other;

        else if ( machNumber <= 1.0 )
@@ -101,11 +106,11 @@ class UpwindMomentumBase 
         if ( machNumber <= -1.0 )

            return density * velocity_main * velocity_other;

        else if ( machNumber <= 0.0 )

            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * velocity_other / speedOfSound;

            return density * speedOfSound / ( 2 * this->gamma ) * ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * velocity_other;

        else if ( machNumber <= 1.0 )

            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( machNumber - 1.0 ) * velocity_other / speedOfSound;

            return density * speedOfSound / ( 2 * this->gamma ) * ( machNumber - 1.0 ) * velocity_other;

        else 

            return 0;

            return 0.0;

      };



      protected:
@@ -118,6 +123,8 @@ class UpwindMomentumBase 
         

         MeshFunctionPointer pressure;



         MeshFunctionPointer density;



};



} //namespace TNL