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

vytvoreny operatory pro jednotlive rovnice 

-nevim jak priradit vytvoreny operator do differentialOperator
-nejsou zatim implementovane operatorove funkce pro fzyikalni veliciny
parent 27dd5b89

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

+51 −0
Original line number Diff line number Diff line
@@ -40,6 +40,57 @@ operator()( const MeshFunction& u, 
    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(); 

      /*

   //rho

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

   typedef typename MeshType::Cell Cell;

   const int size = mesh.template getEntitiesCount< Cell >()/5; 

   return (0.5 * this->tau)( u[ west ] - 2.0 * u[ center ]  + u[ east ] ) 

          - 0.5 * hxInverse * ( u[ west ] * u[ west + 3*size ] - u[ east ] * u[ east + 3*size ] );

   */

   /*

   //rhoVel

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

   typedef typename MeshType::Cell Cell;

   const int size = mesh.template getEntitiesCount< Cell >()/5; 

   return (0.5 * this->tau)( u[ west ] - 2.0 * u[ center ]  + u[ east ] ) 

          - 0.5 * hxInverse * 

          (( u[ west ] * u[ west + 2*size ] + u[ west + 3*size ] ) 

          - (u[ east ] * u[ east + 2*size ] + u[ east + 3*size ] ));

   */

   /*

   //energy

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

   typedef typename MeshType::Cell Cell;

   const int size = mesh.template getEntitiesCount< Cell >()/5; 

   return (0.5 * this->tau)( u[ west ] - 2.0 * u[ center ]  + u[ east ] ) 

          - 0.5 * hxInverse * 

          (( u[ west ] + u[ west + 2*size ] ) * u[ west + size ]  

          - (u[ east ] + u[ east + 2*size ] ) * u[ east + size ] );

   */

   /*

   //vel

   const IndexType& center = entity.getIndex(); 

   typedef typename MeshType::Cell Cell;

   const int size = mesh.template getEntitiesCount< Cell >()/5; 

   return ( fu[ center - 2*size] - fu[ center - 3*size]);

   */

   /*

   //pressure

   const IndexType& center = entity.getIndex(); 

   typedef typename MeshType::Cell Cell;

   const int size = mesh.template getEntitiesCount< Cell >()/5; 

   return ( this->gamma - 1 ) * ( fu[ center - 2*size ] - 0.5 * fu[ center - 3* size ] * fu[ center - 4*size]);

   */



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

   const IndexType& center = entity.getIndex();

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

+4 −0
Original line number Diff line number Diff line
@@ -6,6 +6,10 @@ 
#include <functions/tnlConstantFunction.h>

#include "eulerProblem.h"

#include "LaxFridrichs.h"

#include "LaxFridrichsContinuity.h"

#include "LaxFridrichsMomentum.h"

#include "LaxFridrichsEnergy.h"



#include "eulerRhs.h"

#include "eulerBuildConfigTag.h"

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

+10 −0
Original line number Diff line number Diff line
@@ -25,6 +25,16 @@ class eulerProblem: 
      using typename BaseType::MeshType;

      using typename BaseType::DofVectorType;

      using typename BaseType::MeshDependentDataType;

      typedef tnlMeshFunction<Mesh,Mesh::getMeshDimensions(),RealType>MeshFunction;

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

      tnlSharedVector< RealType, DeviceType, IndexType > rho;

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

+49 −39
Original line number Diff line number Diff line
@@ -106,7 +106,7 @@ setInitialCondition( const tnlParameterContainer& parameters, 
   double eR = ( preR / (gamma - 1) ) + 0.5 * rhoR * velR * velR;

   double x0 = parameters.getParameter< double >( "riemann-border" );

   cout << gamma << " " << rhoL << " " << velL << " " << preL << " " << eL << " " << rhoR << " " << velR << " " << preR << " " << eR << " " << x0 << " " << gamma << endl;

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

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

   this->rho.bind(dofs,0,count);

   this->rhoVel.bind(dofs,count,count);

   this->energy.bind(dofs,2 * count,count);
@@ -248,55 +248,65 @@ getExplicitRHS( const RealType& time, 
    p this->pressure[]

    */

    typedef typename MeshType::Cell Cell;

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

    const RealType& size = 1;//mesh.template getSpaceStepsProducts< -1, 0 >();

    for (long int i = 1; i < count - 1; i++)

       _fu[i] = 1.0 / (2.0*tau) * (this->rho[i-1]+this->rho[i+1]-2.0*this->rho[i])-(1.0/(2.0 * size)) * (this->rho[i+1]

       * this->velocity[i+1] - this->rho[i-1] * this->velocity[i-1]);

       _fu[count-1] = _fu[count-2];

    for (long int i = 1; i < count - 1; i++)

       _fu[count + i] =

       1.0 / (2.0 * tau) * (this->rhoVel[i+1] + this->rhoVel[i+1] - 

       2.0 * this->rhoVel[i])-(1.0 / (2.0 * size)) * ((this->rhoVel[i+1]

       * this->velocity[i + 1] + this->pressure[i + 1]) - (this->rhoVel[i-1] * this->velocity[i - 1] 

       + this->pressure[i - 1]));

       _fu[2*count-1] = _fu[2*count-2];

    for (long int i = 1; i < count - 1; i++)

       _fu[i + 2 * count] = 1.0 / (2.0*tau) * (this->energy[i-1]

       + this->energy[i+1]-2.0*this->energy[i])

       - (1.0/(2.0 * size)) * ((this->energy[i+1]

       + this->pressure[i + 1]) * this->velocity[i + 1] - (this->energy[i-1] + this->pressure[i -1]) * this->velocity[i - 1]);

       _fu[3 * count-1] = _fu[3 * count-2];

    for (long int i = 0; i <= count; i++) //pressure

       this->pressure[i] = (this->gamma - 1 ) * ( this->energy[i] - 0.5 * this->rhoVel[i] * this->velocity[i]);

    for (long int i = 0; i <= count ; i++) //velocity

       this->velocity[i] = this->rhoVel[i]/this->rho[i];

     

   /****

    * If you use an explicit solver like tnlEulerSolver or tnlMersonSolver, you

    * need to implement this method. Compute the right-hand side of

    *

    *   d/dt u(x) = fu( x, u )

    *

    * You may use supporting mesh dependent data if you need.

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

		

	//bind _fu

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

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

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



   MeshFunctionType velocity( mesh, this->velocity );

   MeshFunctionType pressure( mesh, this->pressure );

	//rho

   this->bindDofs( mesh, _u );

   tnlExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;

   MeshFunctionType u( mesh, _u ); 

   MeshFunctionType fu( mesh, _fu ); 

   MeshFunctionType uRho( mesh, _uRho ); 

   MeshFunctionType fuRho( mesh, _fuRho );

   diffrrentialOperatorRho.setTau(tau);

   differentialOperatorRho.setVelocity(velocity) 

   explicitUpdater.template update< typename Mesh::Cell >( time,

                                                           mesh,

                                                           this->differentialOperator,

                                                           this->differentialOperatorRho,

                                                           this->boundaryCondition,

                                                           this->rightHandSide,

                                                           u,

                                                           fu );

                                                           uRho,

                                                           fuRho );

   //rhoVelocity

   MeshFunctionType uRhoVelocity( mesh, _uRhoVelocity ); 

   MeshFunctionType fuRhoVelocity( mesh, _fuRhoVelocity );

   diffrrentialOperatorRhoVelocity.setTau(tau);

   differentialOperatorRhoVelocity.setVelocity(velocity)

   differentialOperatorRhoVelocity.setPressure(pressure) 

   explicitUpdater.template update< typename Mesh::Cell >( time,

                                                           mesh,

                                                           this->differentialOperatorRhoVelocity,

                                                           this->boundaryCondition,

                                                           this->rightHandSide,

                                                           uRhoVelocity,

                                                           fuRhoVelocity );

   

   //energy

   MeshFunctionType uEnergy( mesh, _uEnergy ); 

   MeshFunctionType fuEnergy( mesh, _fuEnergy );

   diffrrentialOperatorEnergy.setTau(tau);

   diffrrentialOperatorEnergy.setPressure(pressure);

   diffrrentialOperatorEnergy.setVelocity(velocity); 

   explicitUpdater.template update< typename Mesh::Cell >( time,

                                                           mesh,

                                                           this->differentialOperatorEnergy,

                                                           this->boundaryCondition,

                                                           this->rightHandSide,

                                                           uEnergy,

                                                           fuEnergy );

   tnlBoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter; 

   boundaryConditionsSetter.template apply< typename Mesh::Cell >( 

      this->boundaryCondition, 

      time + tau, 

       u ); */

       u );

 }



template< typename Mesh,

examples/inviscid-flow/2d/LaxFridrichs_impl.h

+1 −1
Original line number Diff line number Diff line
@@ -95,7 +95,7 @@ operator()( const MeshFunction& u, 






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

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

   const IndexType& center = entity.getIndex(); 

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

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