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

Dodelan merge a 3d verze advekce i eulera, prepsany nazvz vstupn9ch veci,... 

Dodelan merge a 3d verze advekce i eulera, prepsany nazvz vstupn9ch veci, nejsou upravene run skripty
parent a72bda48

examples/advection/LaxFridrichs.h

+21 −1
Original line number Diff line number Diff line
@@ -32,13 +32,18 @@ class LaxFridrichs< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, Real, Index > 
      Real artificalViscosity;

      MeshFunctionType advectionSpeedX;

      MeshFunctionType advectionSpeedY;

      MeshFunctionType advectionSpeedZ;



      void setAdvectionSpeedZ(MeshFunctionType& advectionSpeed)

      {

	   this->advectionSpeedZ.bind(advectionSpeed);

      }



      void setAdvectionSpeedY(MeshFunctionType& advectionSpeed)

      {

	   this->advectionSpeedY.bind(advectionSpeed);

      }





      void setAdvectionSpeedX(MeshFunctionType& advectionSpeed)

      {

	   this->advectionSpeedX.bind(advectionSpeed);
@@ -99,6 +104,14 @@ class LaxFridrichs< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, Real, Index > 
      Real artificalViscosity;

      MeshFunctionType advectionSpeedX;

      MeshFunctionType advectionSpeedY;

      MeshFunctionType advectionSpeedZ;





      void setAdvectionSpeedZ(MeshFunctionType& advectionSpeed)

      {

	   this->advectionSpeedZ.bind(advectionSpeed);

      }





      void setAdvectionSpeedY(MeshFunctionType& advectionSpeed)

      {
@@ -166,6 +179,13 @@ class LaxFridrichs< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Index > 
      Real artificalViscosity;

      MeshFunctionType advectionSpeedX;

      MeshFunctionType advectionSpeedY;

      MeshFunctionType advectionSpeedZ;



      void setAdvectionSpeedZ(MeshFunctionType& advectionSpeed)

      {

	   this->advectionSpeedZ.bind(advectionSpeed);

      }





      void setAdvectionSpeedY(MeshFunctionType& advectionSpeed)

      {

examples/advection/LaxFridrichs_impl.h

+9 −7
Original line number Diff line number Diff line
@@ -267,9 +267,9 @@ operator()( const MeshFunction& u, 
    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 RealType& hxInverse = entity.getMesh().template getSpaceStepsProducts< -2,  0,  0 >(); 

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

   const RealType& hzInverse = 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 >();
@@ -277,9 +277,11 @@ operator()( const MeshFunction& u, 
   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;

   return ( (1.0/6.0) / this->tau ) * this->artificalViscosity * 

          ( u[ west ] + u[ east ] + u[ south ] + u[ north ] + u[ up ] + u [ down ] - 6 * u[ center ] ) -

          (this->advectionSpeedX [ east ] * u[ east ] - this->advectionSpeedX [ west ] * u[ west] ) * hxInverse * 0.5 - 

          (this->advectionSpeedY [ north ] * u[ north ] - this->advectionSpeedY [ south ] * u[ south ] ) * hyInverse * 0.5 -

          (this->advectionSpeedZ [ up ] * u[ up ] - this->advectionSpeedZ [ down ] * u[ down ] ) * hzInverse * 0.5;

}



template< typename MeshReal,
@@ -302,7 +304,7 @@ getLinearSystemRowLength( const MeshType& mesh, 
    * by the Finite difference method.

    */



   //return 2*Dimensions + 1;

   return 2*Dimensions + 1;

}



template< typename MeshReal,

examples/advection/tnlRiemann1DBoundaryConditions.h→examples/advection/Riemann1DBoundaryConditions.h

+39 −28
Original line number Diff line number Diff line 
/*** coppied and changed

//** coppied and changed

/***************************************************************************

                          tnlRiemann1DBoundaryConditions.h  -  description

                             -------------------
@@ -16,21 +16,27 @@ 
 *                                                                         *

 ***************************************************************************/



#ifndef TNLRiemann1DBOUNDARYCONDITIONS_H_

#define TNLRiemann1DBOUNDARYCONDITIONS_H_

#ifndef TNLRIEMANN1DBOUNDARYCONDITIONS_H_

#define TNLRIEMANN1DBOUNDARYCONDITIONS_H_



#include <operators/tnlOperator.h>

#include <functions/tnlConstantFunction.h>

#include <functions/tnlFunctionAdapter.h>

#pragma once



#include <TNL/Operators/Operator.h>

#include <TNL/Functions/Analytic/Constant.h>

#include <TNL/Functions/FunctionAdapter.h>

#include <TNL/Functions/MeshFunction.h>



namespace TNL {

namespace Operators {



template< typename Mesh,

          typename Function = tnlConstantFunction< Mesh::getMeshDimensions(), typename Mesh::RealType >,

          typename Function = Functions::Analytic::Constant< Mesh::getMeshDimensions(), typename Mesh::RealType >,

          int MeshEntitiesDimensions = Mesh::getMeshDimensions(),

          typename Real = typename Mesh::RealType,

          typename Index = typename Mesh::IndexType >

class tnlRiemann1DBoundaryConditions

: public tnlOperator< Mesh,

                      MeshBoundaryDomain,

class Riemann1DBoundaryConditions

: public Operator< Mesh,

                   Functions::MeshBoundaryDomain,

                   MeshEntitiesDimensions,

                   MeshEntitiesDimensions,

                   Real,
@@ -44,21 +50,23 @@ class tnlRiemann1DBoundaryConditions 
      typedef typename MeshType::DeviceType DeviceType;

      typedef Index IndexType;

      

      typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;

      typedef SharedPointer< Mesh > MeshPointer;

      typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;

      typedef typename MeshType::VertexType VertexType;



      static constexpr int getMeshDimensions() { return MeshType::meshDimensions; }



      static void configSetup( tnlConfigDescription& config,

                               const tnlString& prefix = "" )

      static void configSetup( Config::ConfigDescription& config,

                               const String& prefix = "" )

      {

         Function::configSetup( config, prefix );

      }

 

      bool setup( const tnlParameterContainer& parameters,

                  const tnlString& prefix = "" )

      bool setup( const MeshPointer& meshPointer,

                  const Config::ParameterContainer& parameters,

                  const String& prefix = "" )

      {

         return this->function.setup( parameters, prefix );

         return Functions::FunctionAdapter< MeshType, FunctionType >::template setup< MeshPointer >( this->function, meshPointer, parameters, prefix );

      }



      void setFunction( const Function& function )
@@ -117,7 +125,7 @@ class tnlRiemann1DBoundaryConditions 
         typename Matrix::MatrixRow matrixRow = matrix.getRow( entity.getIndex() );

         const IndexType& index = entity.getIndex();

         matrixRow.setElement( 0, index, 1.0 );

         b[ index ] = tnlFunctionAdapter< MeshType, Function >::getValue( this->function, entity, time );

         b[ index ] = Functions::FunctionAdapter< MeshType, Function >::getValue( this->function, entity, time );

      }
@@ -128,12 +136,15 @@ class tnlRiemann1DBoundaryConditions 
   //static_assert( Device::DeviceType == Function::Device::DeviceType );

};





template< typename Mesh,

          typename Function >

ostream& operator << ( ostream& str, const tnlRiemann1DBoundaryConditions< Mesh, Function >& bc )

std::ostream& operator << ( std::ostream& str, const Riemann1DBoundaryConditions< Mesh, Function >& bc )

{

   str << "Riemann1D boundary conditions: vector = " << bc.getVector();

   return str;

}



#endif /* TNLRiemann1DBOUNDARYCONDITIONS_H_ */
 
} // namespace Operators

} // namespace TNL

#endif /* TNLRIEMANN1DBOUNDARYCONDITIONS_H_ */

examples/advection/tnlRiemann2DBoundaryConditions.h→examples/advection/Riemann2DBoundaryConditions.h

+40 −36
Original line number Diff line number Diff line 
/*** coppied and changed

//** coppied and changed

/***************************************************************************

                          tnlRiemann2DBoundaryConditions.h  -  description

                             -------------------
@@ -16,21 +16,27 @@ 
 *                                                                         *

 ***************************************************************************/



#ifndef TNLRiemann2DBOUNDARYCONDITIONS_H_

#define TNLRiemann2DBOUNDARYCONDITIONS_H_

#ifndef TNLRIEMANN2DBOUNDARYCONDITIONS_H_

#define TNLRIEMANN2DBOUNDARYCONDITIONS_H_



#include <operators/tnlOperator.h>

#include <functions/tnlConstantFunction.h>

#include <functions/tnlFunctionAdapter.h>

#pragma once



#include <TNL/Operators/Operator.h>

#include <TNL/Functions/Analytic/Constant.h>

#include <TNL/Functions/FunctionAdapter.h>

#include <TNL/Functions/MeshFunction.h>



namespace TNL {

namespace Operators {



template< typename Mesh,

          typename Function = tnlConstantFunction< Mesh::getMeshDimensions(), typename Mesh::RealType >,

          typename Function = Functions::Analytic::Constant< Mesh::getMeshDimensions(), typename Mesh::RealType >,

          int MeshEntitiesDimensions = Mesh::getMeshDimensions(),

          typename Real = typename Mesh::RealType,

          typename Index = typename Mesh::IndexType >

class tnlRiemann2DBoundaryConditions

: public tnlOperator< Mesh,

                      MeshBoundaryDomain,

class Riemann2DBoundaryConditions

: public Operator< Mesh,

                   Functions::MeshBoundaryDomain,

                   MeshEntitiesDimensions,

                   MeshEntitiesDimensions,

                   Real,
@@ -44,21 +50,23 @@ class tnlRiemann2DBoundaryConditions 
      typedef typename MeshType::DeviceType DeviceType;

      typedef Index IndexType;

      

      typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;

      typedef SharedPointer< Mesh > MeshPointer;

      typedef Containers::Vector< RealType, DeviceType, IndexType> DofVectorType;

      typedef typename MeshType::VertexType VertexType;



      static constexpr int getMeshDimensions() { return MeshType::meshDimensions; }



      static void configSetup( tnlConfigDescription& config,

                               const tnlString& prefix = "" )

      static void configSetup( Config::ConfigDescription& config,

                               const String& prefix = "" )

      {

         Function::configSetup( config, prefix );

      }

 

      bool setup( const tnlParameterContainer& parameters,

                  const tnlString& prefix = "" )

      bool setup( const MeshPointer& meshPointer,

                  const Config::ParameterContainer& parameters,

                  const String& prefix = "" )

      {

         return this->function.setup( parameters, prefix );

         return Functions::FunctionAdapter< MeshType, FunctionType >::template setup< MeshPointer >( this->function, meshPointer, parameters, prefix );

      }



      void setFunction( const Function& function )
@@ -66,11 +74,6 @@ class tnlRiemann2DBoundaryConditions 
         this->function = function;

      }



      void setX0( const RealType& x0 )

      {

         this->x0 = x0;

      }



      Function& getFunction()

      {

         return this->function;
@@ -92,7 +95,7 @@ class tnlRiemann2DBoundaryConditions 
      const auto& neighbourEntities = entity.getNeighbourEntities();

      typedef typename MeshType::Cell Cell;

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

      count = sqrt(count);

      count = std::sqrt(count);

      if( entity.getCoordinates().x() < count * 0.5 )

         return 1;

       else
@@ -125,7 +128,7 @@ class tnlRiemann2DBoundaryConditions 
         typename Matrix::MatrixRow matrixRow = matrix.getRow( entity.getIndex() );

         const IndexType& index = entity.getIndex();

         matrixRow.setElement( 0, index, 1.0 );

         b[ index ] = tnlFunctionAdapter< MeshType, Function >::getValue( this->function, entity, time );

         b[ index ] = Functions::FunctionAdapter< MeshType, Function >::getValue( this->function, entity, time );

      }
@@ -133,17 +136,18 @@ class tnlRiemann2DBoundaryConditions 


      Function function;

 

      RealType x0 = 0.3;

   

   //static_assert( Device::DeviceType == Function::Device::DeviceType );

};





template< typename Mesh,

          typename Function >

ostream& operator << ( ostream& str, const tnlRiemann2DBoundaryConditions< Mesh, Function >& bc )

std::ostream& operator << ( std::ostream& str, const Riemann2DBoundaryConditions< Mesh, Function >& bc )

{

   str << "Riemann2D boundary conditions: vector = " << bc.getVector();

   return str;

}



#endif /* TNLRiemann2DBOUNDARYCONDITIONS_H_ */
 
} // namespace Operators

} // namespace TNL

#endif /* RIEMANN2DBOUNDARYCONDITIONS_H_ */

examples/advection/advection.h

+15 −32
Original line number Diff line number Diff line
@@ -8,8 +8,8 @@ 
#include "LaxFridrichs.h"

#include "advectionRhs.h"

#include "advectionBuildConfigTag.h"

#include "tnlRiemann1DBoundaryConditions.h"

#include "tnlRiemann2DBoundaryConditions.h"

#include "Riemann1DBoundaryConditions.h"

#include "Riemann2DBoundaryConditions.h"



using namespace TNL;
@@ -31,37 +31,28 @@ template< typename ConfigTag >class advectionConfig 
      static void configSetup( Config::ConfigDescription & config )

      {

         config.addDelimiter( "advection settings:" );

<<<<<<< HEAD

         config.addEntry< tnlString >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet");

            config.addEntryEnum< tnlString >( "dirichlet" );

            config.addEntryEnum< tnlString >( "neumann" );

            config.addEntryEnum< tnlString >( "riemann1D" );

            config.addEntryEnum< tnlString >( "riemann2D" );

         config.addEntry< double >( "boundary-conditions-constant", "This sets a value in case of the constant boundary conditions." );

	 config.addEntry< double >( "artifical-viscosity", "This sets value of artifical viscosity (default 1)", 1.0);

	 config.addEntry< tnlString >( "begin", "choose begin type", "sin");

	    config.addEntryEnum< tnlString >( "exp");

	    config.addEntryEnum< tnlString >( "exp_square");

	    config.addEntryEnum< tnlString >( "square");

	    config.addEntryEnum< tnlString >( "riemann");

=======

         config.addEntry< String >( "boundary-conditions-type", "Choose the boundary conditions type.", "dirichlet");

            config.addEntryEnum< String >( "dirichlet" );

            config.addEntryEnum< String >( "neumann" );

            config.addEntryEnum< String >( "riemann1D" );

            config.addEntryEnum< String >( "riemann2D" );

         config.addEntry< double >( "boundary-conditions-constant", "This sets a value in case of the constant boundary conditions." );

	 config.addEntry< double >( "artifical-viscosity", "This sets value of artifical viscosity (default 1)", 1.0);

	 config.addEntry< String >( "begin", "choose begin type", "sin");

	    config.addEntryEnum< String >( "sin");

	    config.addEntryEnum< String >( "sin_square");

>>>>>>> develop

	    config.addEntryEnum< String >( "exp");

	    config.addEntryEnum< String >( "exp_square");

	    config.addEntryEnum< String >( "square");

	    config.addEntryEnum< String >( "riemann");

	 config.addEntry< double >( "advection-speedX", "This sets value of advection speed in X direction (default 1)" , 1.0);

	 config.addEntry< double >( "advection-speedY", "This sets value of advection speed in Y direction (default 1)" , 1.0);

	 config.addEntry< double >( "advection-speedZ", "This sets value of advection speed in Z direction (default 1)" , 1.0);

	 config.addEntry< String >( "move", "choose movement type", "advection");

	    config.addEntryEnum< String >( "advection");

	    config.addEntryEnum< String >( "rotation");

	 config.addEntry< int >( "dimension", "choose movement typeproblem dimension", 1);

	    config.addEntryEnum< int >( 1 );

	    config.addEntryEnum< int >( 2 );

            config.addEntryEnum< int >( 3 );

	 config.addEntry< double >( "realSize", "Real size of scheme", 1.0);



         /****
@@ -121,36 +112,28 @@ class advectionSetter 
             SolverStarter solverStarter;

             return solverStarter.template run< Problem >( parameters );

          }

<<<<<<< HEAD

          if( boundaryConditionsType == "riemann1D" )

          {

             typedef tnlRiemann1DBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;

             typedef Operators::Riemann1DBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;

             typedef advectionProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;

             SolverStarter solverStarter;

             return solverStarter.template run< Problem >( parameters );

          }

          if( boundaryConditionsType == "riemann2D" )

          {

             typedef tnlRiemann2DBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;

             typedef Operators::Riemann2DBoundaryConditions< MeshType, MeshFunction, MeshType::getMeshDimensions(), Real, Index > BoundaryConditions;

             typedef advectionProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;

             SolverStarter solverStarter;

             return solverStarter.template run< Problem >( parameters );

          }

          if( boundaryConditionsType == "neumann" )

          {

             typedef tnlNeumannBoundaryConditions< MeshType, MeshFunction, Real, Index > BoundaryConditions;

             typedef advectionProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;

             SolverStarter solverStarter;

             return solverStarter.template run< Problem >( parameters );

          }

=======

             typedef Operators::NeumannBoundaryConditions< MeshType, MeshFunction, Real, Index > BoundaryConditions;

             typedef advectionProblem< MeshType, BoundaryConditions, RightHandSide, ApproximateOperator > Problem;

             SolverStarter solverStarter;

             return solverStarter.template run< Problem >( parameters );

>>>>>>> develop

          }



      return true;}

};



int main( int argc, char* argv[] )