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

vzresena 1D hustota

parent 888f2a4d
Loading
Loading
Loading
Loading
+14 −9
Original line number Diff line number Diff line
@@ -65,7 +65,7 @@ class UpwindContinuityBase
          this->velocity = velocity;
      };

      const RealType& positiveDensityFlux( const RealType& density, const RealType& velocity, const RealType& pressure )
      RealType positiveDensityFlux( const RealType& density, const RealType& velocity, const RealType& pressure ) const
      {
         const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
         const RealType& machNumber = velocity / speedOfSound;
@@ -77,9 +77,9 @@ class UpwindContinuityBase
            return density * speedOfSound / ( 2 * this->gamma ) * ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 );
        else 
            return density * velocity;
      }
      };

      const RealType& negativeDensityFlux( const RealType& density, const RealType& velocity, const RealType& pressure )
      RealType negativeDensityFlux( const RealType& density, const RealType& velocity, const RealType& pressure ) const
      {
         const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
         const RealType& machNumber = velocity / speedOfSound;
@@ -91,7 +91,12 @@ class UpwindContinuityBase
            return density * speedOfSound / ( 2 * this->gamma ) * ( machNumber - 1.0 );
        else 
            return 0.0;
      }
      };
      
      RealType multiply (const RealType& a, const RealType& b ) const
      {
         return a * b;
      };
      

      protected:
@@ -162,11 +167,11 @@ class UpwindContinuity< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, In
         const RealType& velocity_x_east   = this->velocity.template getData< DeviceType >()[ 0 ].template getData< DeviceType >()[ east ];

         return -hxInverse * (
                                   UpwindContinuity::positiveDensityFlux( u[ center ], velocity_x_center, pressure_center )
                                -  UpwindContinuity::positiveDensityFlux( u[ west   ], velocity_x_west  , pressure_west   )
                                -  UpwindContinuity::negativeDensityFlux( u[ center ], velocity_x_center, pressure_center )
                                +  UpwindContinuity::negativeDensityFlux( u[ east   ], velocity_x_east  , pressure_east   )
//                                  u[ center ] / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_x_center + std::sqrt( this->gamma * pressure_center / u[ center ] ) )
                                   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   )
//                                  u[ center ] / ( 2 * this->gamma ) * ( ( 2 * this->gamma - 1 ) * velocity_x_center + std::sqrt( this->multiply( 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   ] ) )
+15 −15
Original line number Diff line number Diff line
@@ -67,7 +67,7 @@ class UpwindEnergyBase
      void setArtificialViscosity( const RealType& artificialViscosity )
      {
         this->artificialViscosity = artificialViscosity;
      }
      };

      const RealType& positiveEnergyFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other1, const RealType& velocity_other2, const RealType& pressure )
      {
@@ -78,44 +78,44 @@ class UpwindEnergyBase
        else if ( machNumber <= 0.0 )
            return density * speedOfSound / ( 2 * this->gamma ) 
                 * ( 
                     ( speedOfSound * ( machNumber + 1.0 ) * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) );
                   + ( speedOfSound * speedOfSound * machNumber * ( - 1.0 - machNumber )
                     ( speedOfSound * ( machNumber + 1.0 ) * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
                   + ( speedOfSound * speedOfSound * machNumber * ( - 1.0 - machNumber ) )
                   + ( speedOfSound * speedOfSound * ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
                   );
        else if ( machNumber <= 1.0 )
            return density * speedOfSound / ( 2 * this->gamma ) 
                 * ( 
                     ( speedOfSound * ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) ) * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) );
                   + ( speedOfSound * speedOfSound * machNumber * ( machNumber + 1.0 )
                     ( speedOfSound * ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
                   + ( speedOfSound * speedOfSound * machNumber * ( machNumber + 1.0 ) )
                   + ( speedOfSound * speedOfSound * ( machNumber + 1.0 ) / ( this->gamma - 1.0 ) )
                   );
        else   
            return ( pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 );
      }
            return ( pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) );
      };

      const RealType& negativeEnergyFlux( const RealType& density, const RealType& velocity, const RealType& pressure )
      const RealType& negativeEnergyFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other1, const RealType& velocity_other2, const RealType& pressure )
      {
         const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
         const RealType& machNumber = velocity / speedOfSound;
         if ( machNumber <= -1.0 )
            return ( pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 );
            return ( pressure / ( this->gamma - 1.0 ) + 0.5 * density * ( velocity_main * velocity_main + velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) );
        else if ( machNumber <= 0.0 )
            return density * speedOfSound / ( 2 * this->gamma ) 
                 * ( 
                     ( speedOfSound * ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 )  * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) );
                   + ( speedOfSound * speedOfSound * machNumber * ( machNumber - 1.0 )
                     ( speedOfSound * ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 )  * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
                   + ( speedOfSound * speedOfSound * machNumber * ( machNumber - 1.0 ) )
                   + ( speedOfSound * speedOfSound * ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
                   )
                   );
        else if ( machNumber <= 1.0 )
            return density * speedOfSound / ( 2 * this->gamma ) 
                 * ( 
                     ( speedOfSound * ( machNumber - 1.0 ) * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) );
                   + ( speedOfSound * speedOfSound * machNumber * ( 1.0 - machNumber )
                     ( speedOfSound * ( machNumber - 1.0 ) * speedOfSound * ( machNumber * machNumber + ( velocity_other1 * velocity_other1 + velocity_other2 * velocity_other2 ) / ( speedOfSound * speedOfSound ) ) )
                   + ( speedOfSound * speedOfSound * machNumber * ( 1.0 - machNumber ) )
                   + ( speedOfSound * speedOfSound * ( machNumber - 1.0 ) / ( this->gamma - 1.0 ) )
                   );
        else 
            return 0.0;
      }      
      };      

      protected:
         
+12 −12
Original line number Diff line number Diff line
@@ -59,26 +59,26 @@ class UpwindMomentumBase
         if ( machNumber <= -1.0 )
            return 0;
        else if ( machNumber <= 0.0 )
            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( machNumber + 1.0 ) * machNumber + ( - machNumber - 1.0 );
            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( machNumber + 1.0 ) * machNumber + ( - machNumber - 1.0 ) );
        else if ( machNumber <= 1.0 )
            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) * machNumber + (machNumber + 1.0 );
            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) * machNumber + (machNumber + 1.0 ) );
        else 
            return density * velocity * velocity + pressure;
      }
      };

      const RealType& negativeMainMomentumFlux( const RealType& density, const RealType& velocity, const RealType& pressure )
      {
         const RealType& speedOfSound = std::sqrt( this->gamma * pressure / density );
         const RealType& machNumber = velocity / speedOfSound;
         if ( machNumber <= -1.0 )
            return density * velocity;
            return density * velocity * velocity + pressure;
        else if ( machNumber <= 0.0 )
            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * machNumber + (machNumber - 1.0 );
            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * machNumber + (machNumber - 1.0 ) );
        else if ( machNumber <= 1.0 )
            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( machNumber - 1.0 ) * machNumber + ( - machNumber + 1.0 );
            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( machNumber - 1.0 ) * machNumber + ( - machNumber + 1.0 ) );
        else 
            return density * velocity * velocity + pressure;
      }
            return 0; 
      };

      const RealType& positiveOtherMomentumFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other, const RealType& pressure )
      {
@@ -92,7 +92,7 @@ class UpwindMomentumBase
            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( ( 2.0 * this->gamma - 1.0 ) * machNumber + 1.0 ) * velocity_other / speedOfSound;
        else 
            return density * velocity_main * velocity_other;
      }
      };

      const RealType& negativeOtherMomentumFlux( const RealType& density, const RealType& velocity_main, const RealType& velocity_other, const RealType& pressure )
      {
@@ -101,12 +101,12 @@ 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 * speedOfSound / ( 2 * this->gamma ) * ( ( 2.0 * this->gamma - 1.0 ) * machNumber - 1.0 ) * velocity_other / speedOfSound;
        else if ( machNumber <= 1.0 )
            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( machNumber + 1.0 ) * velocity_other / speedOfSound;
            return density * speedOfSound * speedOfSound / ( 2 * this->gamma ) * ( machNumber - 1.0 ) * velocity_other / speedOfSound;
        else 
            return 0;
      }
      };

      protected: