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

Adding vector field. 

Refactoring the Lax-Fridrichs scheme for the transport equation.
parent c6344062

examples/advection/LaxFridrichs.h

+166 −113
Original line number Diff line number Diff line
@@ -3,12 +3,14 @@ 


#include <TNL/Containers/Vector.h>

#include <TNL/Meshes/Grid.h>

#include <TNL/Functions/VectorField.h>



namespace TNL {



template< typename Mesh,

          typename Real = typename Mesh::RealType,

          typename Index = typename Mesh::IndexType >

          typename Index = typename Mesh::IndexType,

          typename VelocityFunction = Functions::MeshFunction< Mesh > >

class LaxFridrichs

{

};
@@ -17,36 +19,35 @@ template< typename MeshReal, 
          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

class LaxFridrichs< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, Real, Index >

          typename Index,

          typename VelocityFunction >

class LaxFridrichs< Meshes::Grid< 1, MeshReal, Device, MeshIndex >, Real, Index, VelocityFunction >

{

   public:

      

      typedef Meshes::Grid< 1, MeshReal, Device, MeshIndex > MeshType;

      static const int Dimensions = MeshType::getMeshDimensions();

      typedef typename MeshType::CoordinatesType CoordinatesType;

      typedef Real RealType;

      typedef Device DeviceType;

      typedef Index IndexType;

      typedef Functions::MeshFunction< MeshType > MeshFunctionType;

      enum { Dimensions = MeshType::getMeshDimensions() };

      Real tau;

      Real artificalViscosity;

      Real advectionSpeedX;

      Real advectionSpeedY;



      void setAdvectionSpeedY(const Real& advectionSpeed)

      {

	   this->advectionSpeedY = advectionSpeed;

      }

      typedef VelocityFunction VelocityFunctionType;

      typedef Functions::VectorField< Dimensions, VelocityFunctionType > VelocityFieldType;

      

      LaxFridrichs() : artificialViscosity( 1.0 ) {}

      

      void setAdvectionSpeedX(const Real& advectionSpeed)

      bool setup( const Config::ParameterContainer& parameters,

                  const String& prefix = "" )

      {

	   this->advectionSpeedX = advectionSpeed;

         return true;

      }



      static String getType();

      

      void setViscosity(const Real& artificalViscosity)

      {

	   this->artificalViscosity = artificalViscosity;

         this->artificialViscosity = artificalViscosity;

      }

      

      void setTau(const Real& tau)
@@ -54,66 +55,78 @@ class LaxFridrichs< Meshes::Grid< 1,MeshReal, Device, MeshIndex >, Real, Index > 
          this->tau = tau;

      };



      static String getType();

      VelocityFieldType& getVelocityField()

      {

         return this->velocityField;

      }

      

      const VelocityFieldType& getVelocityField() const

      {

         return this->velocityField;

      }

      

      template< typename MeshFunction, typename MeshEntity >

      __cuda_callable__

      Real operator()( const MeshFunction& u,

                       const MeshEntity& entity,

                       const RealType& time = 0.0 ) const;

                       const RealType& time = 0.0 ) const

      {

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



         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 Functions::FunctionAdapter< MeshType, VelocityFunctionType > FunctionAdapter;

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

                FunctionAdapter::getValue( this->velocityField[ 0 ], entity, time ) * ( u[ east ] - u[west] ) * hxInverse * 0.5;

      }

      

      template< typename MeshEntity >

      __cuda_callable__

      Index getLinearSystemRowLength( const MeshType& mesh,

                                      const IndexType& index,

                                      const MeshEntity& entity ) const;

   protected:

            

      template< typename MeshEntity, typename Vector, typename MatrixRow >

      __cuda_callable__

      void updateLinearSystem( const RealType& time,

                               const RealType& tau,

                               const MeshType& mesh,

                               const IndexType& index,

                               const MeshEntity& entity,

                               const MeshFunctionType& u,

                               Vector& b,

                               MatrixRow& matrixRow ) const;

      RealType tau;

      

      RealType artificialViscosity;

      

      VelocityFieldType velocityField;

};



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

class LaxFridrichs< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, Real, Index >

          typename Index,

          typename VelocityFunction >

class LaxFridrichs< Meshes::Grid< 2, MeshReal, Device, MeshIndex >, Real, Index, VelocityFunction >

{

   public:

      

      typedef Meshes::Grid< 2, MeshReal, Device, MeshIndex > MeshType;

      static const int Dimensions = MeshType::getMeshDimensions();

      typedef typename MeshType::CoordinatesType CoordinatesType;

      typedef Real RealType;

      typedef Device DeviceType;

      typedef Index IndexType;

      typedef Functions::MeshFunction< MeshType > MeshFunctionType;

      enum { Dimensions = MeshType::getMeshDimensions() };

      Real tau;

      Real artificalViscosity;

      Real advectionSpeedX;

      Real advectionSpeedY;



      void setAdvectionSpeedY(const Real& advectionSpeed)

      {

	   this->advectionSpeedY = advectionSpeed;

      }

      typedef VelocityFunction VelocityFunctionType;

      typedef Functions::VectorField< Dimensions, VelocityFunctionType > VelocityFieldType;

      

      LaxFridrichs()

         : artificialViscosity( 1.0 ) {}      

      

      void setAdvectionSpeedX(const Real& advectionSpeed)

      bool setup( const Config::ParameterContainer& parameters,

                  const String& prefix = "" )

      {

	   this->advectionSpeedX = advectionSpeed;

         return true;

      }



      static String getType();

      

      void setViscosity(const Real& artificalViscosity)

      {

	   this->artificalViscosity = artificalViscosity;

         this->artificialViscosity = artificalViscosity;

      }

      

      void setTau(const Real& tau)
@@ -121,66 +134,84 @@ class LaxFridrichs< Meshes::Grid< 2,MeshReal, Device, MeshIndex >, Real, Index > 
          this->tau = tau;

      };



      static String getType();

      VelocityFieldType& getVelocityField()

      {

         return this->velocityField;

      }

      

      const VelocityFieldType& getVelocityField() const

      {

         return this->velocityField;

      }

      

      template< typename MeshFunction, typename MeshEntity >

      __cuda_callable__

      Real operator()( const MeshFunction& u,

                       const MeshEntity& entity,

                       const RealType& time = 0.0 ) const;

                       const RealType& time = 0.0 ) const

      {

         static_assert( MeshEntity::entityDimensions == 2, "Wrong mesh entity dimensions." ); 

         static_assert( MeshFunction::getEntitiesDimensions() == 2, "Wrong preimage function" ); 

         const typename MeshEntity::template NeighbourEntities< 2 >& neighbourEntities = entity.getNeighbourEntities(); 



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

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

         

         const IndexType& center = entity.getIndex();

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

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

         const IndexType& north = neighbourEntities.template getEntityIndex<  0,  1 >(); 

         const IndexType& south = neighbourEntities.template getEntityIndex<  0, -1 >(); 

         

         typedef Functions::FunctionAdapter< MeshType, VelocityFunctionType > FunctionAdapter;

         return ( 0.25 / this->tau ) * this->artificialViscosity * ( u[ west ] + u[ east ] + u[ north ] + u[ south ] - 4.0 * u[ center ] ) -

                0.5 * ( FunctionAdapter::getValue( this->velocityField[ 0 ], entity, time ) * ( u[ east ] - u[ west ] ) * hxInverse +

                        FunctionAdapter::getValue( this->velocityField[ 1 ], entity, time ) * ( u[ north ] - u[ south ] ) * hyInverse );         

      }

      

      template< typename MeshEntity >

      __cuda_callable__

      Index getLinearSystemRowLength( const MeshType& mesh,

                                      const IndexType& index,

                                      const MeshEntity& entity ) const;

   protected:

            

      template< typename MeshEntity, typename Vector, typename MatrixRow >

      __cuda_callable__

      void updateLinearSystem( const RealType& time,

                               const RealType& tau,

                               const MeshType& mesh,

                               const IndexType& index,

                               const MeshEntity& entity,

                               const MeshFunctionType& u,

                               Vector& b,

                               MatrixRow& matrixRow ) const;

      RealType tau;

      

      RealType artificialViscosity;

      

      VelocityFieldType velocityField;

};



template< typename MeshReal,

          typename Device,

          typename MeshIndex,

          typename Real,

          typename Index >

class LaxFridrichs< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Index >

          typename Index,

          typename VelocityFunction >

class LaxFridrichs< Meshes::Grid< 3, MeshReal, Device, MeshIndex >, Real, Index, VelocityFunction >

{

   public:

      

      typedef Meshes::Grid< 3, MeshReal, Device, MeshIndex > MeshType;

      static const int Dimensions = MeshType::getMeshDimensions();

      typedef typename MeshType::CoordinatesType CoordinatesType;

      typedef Real RealType;

      typedef Device DeviceType;

      typedef Index IndexType;

      typedef Functions::MeshFunction< MeshType > MeshFunctionType;

      enum { Dimensions = MeshType::getMeshDimensions() };

      Real tau;

      Real artificalViscosity;

      Real advectionSpeedX;

      Real advectionSpeedY;

      typedef VelocityFunction VelocityFunctionType;

      typedef Functions::VectorField< Dimensions, VelocityFunctionType > VelocityFieldType;

      

      void setAdvectionSpeedY(const Real& advectionSpeed)

      {

	   this->advectionSpeedY = advectionSpeed;

      }

      LaxFridrichs()

         : artificialViscosity( 1.0 ) {}      

      



      void setAdvectionSpeedX(const Real& advectionSpeed)

      bool setup( const Config::ParameterContainer& parameters,

                  const String& prefix = "" )

      {

	   this->advectionSpeedX = advectionSpeed;

         return true;

      }



      static String getType();

      

      void setViscosity(const Real& artificalViscosity)

      {

	   this->artificalViscosity = artificalViscosity;

         this->artificialViscosity = artificalViscosity;

      }

      

      void setTau(const Real& tau)
@@ -188,35 +219,57 @@ class LaxFridrichs< Meshes::Grid< 3,MeshReal, Device, MeshIndex >, Real, Index > 
          this->tau = tau;

      };



      static String getType();

      VelocityFieldType& getVelocityField()

      {

         return this->velocityField;

      }

      

      const VelocityFieldType& getVelocityField() const

      {

         return this->velocityField;

      }

      

      template< typename MeshFunction, typename MeshEntity >

      __cuda_callable__

      Real operator()( const MeshFunction& u,

                       const MeshEntity& entity,

                       const RealType& time = 0.0 ) const;

                       const RealType& time = 0.0 ) const

      {

         static_assert( MeshEntity::entityDimensions == 3, "Wrong mesh entity dimensions." ); 

         static_assert( MeshFunction::getEntitiesDimensions() == 3, "Wrong preimage function" ); 

         const typename MeshEntity::template NeighbourEntities< 3 >& neighbourEntities = entity.getNeighbourEntities(); 



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

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

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

         const IndexType& center = entity.getIndex();

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

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

         const IndexType& north = neighbourEntities.template getEntityIndex<  0,  1,  0 >(); 

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

         

         typedef Functions::FunctionAdapter< MeshType, VelocityFunctionType > FunctionAdapter;

         return ( 0.25 / this->tau ) * this->artificialViscosity * ( u[ west ] + u[ east ] + u[ north ] + u[ south ] + u[ up ] + u[ down ] - 6.0 * u[ center ] ) -

                0.5 * ( FunctionAdapter::getValue( this->velocityField[ 0 ], entity, time ) * ( u[ east ] - u[ west ] ) * hxInverse +

                        FunctionAdapter::getValue( this->velocityField[ 1 ], entity, time ) * ( u[ north ] - u[ south ] ) * hyInverse +

                        FunctionAdapter::getValue( this->velocityField[ 2 ], entity, time ) * ( u[ up ] - u[ down ] ) * hzInverse );         

      }

      

      template< typename MeshEntity >

      __cuda_callable__

      Index getLinearSystemRowLength( const MeshType& mesh,

                                      const IndexType& index,

                                      const MeshEntity& entity ) const;

   protected:

            

      template< typename MeshEntity, typename Vector, typename MatrixRow >

      __cuda_callable__

      void updateLinearSystem( const RealType& time,

                               const RealType& tau,

                               const MeshType& mesh,

                               const IndexType& index,

                               const MeshEntity& entity,

                               const MeshFunctionType& u,

                               Vector& b,

                               MatrixRow& matrixRow ) const;

      RealType tau;

      

      RealType artificialViscosity;

      

      VelocityFieldType velocityField;

};





} // namespace TNL





#include "LaxFridrichs_impl.h"

//#include "LaxFridrichs_impl.h"



#endif	/* LaxFridrichs_H */

examples/advection/advectionProblem_impl.h

+16 −100
Original line number Diff line number Diff line
@@ -55,6 +55,14 @@ setup( const MeshPointer& meshPointer, 
       const Config::ParameterContainer& parameters,

       const String& prefix )

{

   this->velocityType = parameters.getParameter< String >( "move" );

   const double artificalViscosity = parameters.getParameter< double >( "artifical-viscosity" );

   differentialOperatorPointer->setViscosity(artificalViscosity);

   const double advectionSpeedX = parameters.getParameter< double >( "advection-speedX" );

   differentialOperatorPointer->setAdvectionSpeedX(advectionSpeedX);

   const double advectionSpeedY = parameters.getParameter< double >( "advection-speedY" );

   differentialOperatorPointer->setAdvectionSpeedY(advectionSpeedY);



   if( ! this->boundaryConditionPointer->setup( meshPointer, parameters, prefix + "boundary-conditions-" ) ||

       ! this->rightHandSidePointer->setup( parameters, prefix + "right-hand-side-" ) )

      return false;
@@ -82,9 +90,11 @@ template< typename Mesh, 
          typename DifferentialOperator >

void

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::

bindDofs( const MeshPointer& mesh,

bindDofs( const MeshPointer& meshPointer,

          DofVectorPointer& dofVector )

{

   const IndexType dofs = meshPointer->template getEntitiesCount< typename MeshType::Cell >();

   this->uPointer->bind( meshPointer, dofVector );

}



template< typename Mesh,
@@ -94,87 +104,17 @@ template< typename Mesh, 
bool

advectionProblem< Mesh, BoundaryCondition, RightHandSide, DifferentialOperator >::

setInitialCondition( const Config::ParameterContainer& parameters,

                     const MeshPointer& mesh,

                     const MeshPointer& meshPointer,

                     DofVectorPointer& dofs,

                     MeshDependentDataPointer& meshDependentData )

{

   std::cout << "vaules adding";

   typedef typename MeshType::Cell Cell;

   int dimensions = parameters.getParameter< int >( "dimension" );

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

   const RealType& size = parameters.getParameter< double >( "realSize" ) / ::pow(count, 1.0/dimensions);

   const String& beginChoice = parameters.getParameter< String >( "begin" );

   std::cout << beginChoice << " " << dimensions << "   " << size << "   " << count << "   "<< 1/dimensions << std::endl;

   //getchar();

   /*if (beginChoice == "step")

   {

	   if (dimensions == 1)

	   {

         for( IndexType i = 0; i < count; i++)

         {

            if( i > 0.45*count  && i<0.55*count )

               dofs->setElement( i, 1.0 );

            else

               dofs->setElement( i, 0.0 );

         }

		}

      else if (dimensions == 2)

	       {

                   count = std::sqrt(count);

		   double expValue;

		   for (IndexType i = 0; i < count-1; i++)

                      for (IndexType j = 0; j < count-1; j++)

		      {

			expValue = std::exp(-std::pow(size*i-2,2)-std::pow(size*j-2,2));

			if( (i>0.4*count) && (i<0.5*count) && (j>0.4*count) && (j<0.5*count) )

                           constantFunction=1;

                        else constantFunction=0;

			if( expValue>constantFunction)

                            ( *dofs )[i * count + j] = expValue;

                        else ( *dofs )[i * count + j] = constantFunction;

		      };

		};

       }

   else if (beginChoice == "sin")

      {

	   if (dimensions == 1)

	      {

		   ( *dofs )[0] = 0;

		   for (IndexType i = 1; i < count-2; i++)

		   {

			( *dofs )[i] = std::exp(-std::pow(size*i-2,2));

		   };

		   ( *dofs )[count-1] = 0;

		}

	    else if (dimensions == 2)

	       {

                   count = ::sqrt(count);

		   for (IndexType i = 1; i < count-1; i++)

		      for (IndexType j = 1; j < count-1; j++)

		      {

			   ( *dofs )[i * count + j] = std::exp(-std::pow(size*i-2,2)-std::pow(size*j-2,2));

		      };

		};

     };*/

   //setting velocity field

   //std::cout << *dofs << std::endl;

   //getchar();

   /*const String& initialConditionFile = parameters.getParameter< String >( "initial-condition" );

   if( ! dofs.load( initialConditionFile ) )

   this->bindDofs( meshPointer, dofs );

   const String& initialConditionFile = parameters.getParameter< String >( "initial-condition" );

   if( ! this->uPointer->boundLoad( initialConditionFile ) )

   {

      std::cerr << "I am not able to load the initial condition from the file " << initialConditionFile << "." << std::endl;

      return false;

   }

   return true;*/

   dofs->save( "dofs.tnl" );

   this->velocityType = parameters.getParameter< String >( "move" );

   const double artificalViscosity = parameters.getParameter< double >( "artifical-viscosity" );

   differentialOperatorPointer->setViscosity(artificalViscosity);

   const double advectionSpeedX = parameters.getParameter< double >( "advection-speedX" );

   differentialOperatorPointer->setAdvectionSpeedX(advectionSpeedX);

   const double advectionSpeedY = parameters.getParameter< double >( "advection-speedY" );

   differentialOperatorPointer->setAdvectionSpeedY(advectionSpeedY);

   std::cout << "vaules added";

   return true;

}
@@ -250,29 +190,6 @@ getExplicitRHS( const RealType& time, 
    */

   typedef typename MeshType::Cell Cell;

   int count = ::sqrt(mesh->template getEntitiesCount< Cell >());

//   const RealType& size = parameters.getParameter< double >( "realSize" ) / ::pow(count, 0.5);

/*   if (this->velocityType == "rotation")

   {

      double radius;

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

         for (int j =1; j < count; j++)

            {

               radius = ::sqrt(pow(i-1-(count/2.0),2) + ::pow(j-1-(count/2.0),2));

            if (radius != 0.0)

               _fu[(i-1)*count+j-1] =(0.25*tau)*differentialOperator.artificalViscosity*			//smoothening part

               (_u[(i-1)*count-2+j]+_u[(i-1)*count+j]+

               _u[i*count+j-1]+_u[(i-2)*count+j-1]- 

               4.0*_u[(i-1)*count+j-1])

               -((1.0/(2.0*count))*differentialOperator.advectionSpeedX						//X addition

               *radius*(-1)*((j-1-(count/2.0))/radius)

	       *(_u[(i-1)*count+j]-_u[(i-1)*count+j-2])) 

	       -((1.0/(2.0*count))*differentialOperator.advectionSpeedY						//Y addition

               *radius*((i-1-(count/2.0))/radius)

	       *(_u[i*count+j-1]-_u[(i-2)*count+j-1]))

            ;}

  }

   else if (this->velocityType == "advection")

*/  { 

   this->bindDofs( mesh, _u );

   Solvers::PDE::ExplicitUpdater< Mesh, MeshFunctionType, DifferentialOperator, BoundaryCondition, RightHandSide > explicitUpdater;

   SharedPointer< MeshFunctionType > u( mesh, _u );
@@ -291,7 +208,6 @@ getExplicitRHS( const RealType& time, 
      time + tau, 

       u ); */

}

}

template< typename Mesh,

          typename BoundaryCondition,

          typename RightHandSide,

examples/advection/tnl-run-advection

+23 −14
Original line number Diff line number Diff line 
#!/usr/bin/env bash



tnl-grid-setup --dimensions 1 \

tnl-grid-setup --dimensions 2 \

               --origin-x 0.0 \

               --proportions-x 10.0 \

               --size-x 500 \

               --origin-x 0.0 \

               --proportions-x 1.0 \

               --proportions-y 1.0 \

               --size-x 100 \

               --size-y 100 

 

tnl-init --test-function heaviside-of-paraboloid \

         --output-file init.tnl \

         --x-center 0.5 \

         --y-center 0.5 \

         --z-center 0.5 \

         --radius 0.1



#tnl-init --test-function sin-wave \

#         --output-file init.tnl

tnl-advection-dbg --time-discretisation explicit \

	      --time-step 1.0e-3 \

              --boundary-conditions-constant 0 \

              --discrete-solver merson \

tnl-advection-dbg --initial-condition init.tnl \

                  --time-discretisation explicit \

                  --time-step 1.0e-5 \

                  --boundary-conditions-constant 1.0 \

                  --discrete-solver euler \

                  --snapshot-period 0.1 \

                  --final-time 1.0 \

                  --artifical-viscosity 1.0 \

	      --begin sin \

	          --begin sin_square \

	          --advection-speedX 2.0 \



tnl-view --mesh mesh.tnl --input-files *tnl

src/TNL/Functions/VectorField.h

0 → 100644
+96 −0

File added.

Preview size limit exceeded, changes collapsed.

src/TNL/Operators/Analytic/Shift.h

+2 −2
Original line number Diff line number Diff line
@@ -20,8 +20,8 @@ namespace Analytic { 
   

   

template< typename Function >

class Shift : public Functions::Domain< Function::getDomainDimenions(), 

                                        Function::getDomainTyep() >

class Shift : public Functions::Domain< Function::getDomainDimensions(), 

                                        Function::getDomainType() >

{

   public: