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

royjeta advekce, u prirazeni meshfunction chzba viz email.

parent 40e37b4a

examples/advection/LaxFridrichs.h

+18 −18
Original line number Diff line number Diff line
@@ -27,22 +27,22 @@ class LaxFridrichs< tnlGrid< 1,MeshReal, Device, MeshIndex >, Real, Index > 
      typedef tnlMeshFunction< MeshType > MeshFunctionType;

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

      Real tau;

      double artificalViscosity;

      double advectionSpeedX;

      double advectionSpeedY;

      Real artificalViscosity;

      Real advectionSpeedX;

      Real advectionSpeedY;



      void setAdvectionSpeedY(const double advectionSpeed)

      void setAdvectionSpeedY(const Real& advectionSpeed)

      {

	   this->advectionSpeedY = advectionSpeed;

      }





      void setAdvectionSpeedX(const double advectionSpeed)

      void setAdvectionSpeedX(const Real& advectionSpeed)

      {

	   this->advectionSpeedX = advectionSpeed;

      }



      void setViscosity(const double artificalViscosity)

      void setViscosity(const Real& artificalViscosity)

      {

	   this->artificalViscosity = artificalViscosity;

      }
@@ -94,22 +94,22 @@ class LaxFridrichs< tnlGrid< 2,MeshReal, Device, MeshIndex >, Real, Index > 
      typedef tnlMeshFunction< MeshType > MeshFunctionType;

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

      Real tau;

      double artificalViscosity;

      double advectionSpeedX;

      double advectionSpeedY;

      Real artificalViscosity;

      Real advectionSpeedX;

      Real advectionSpeedY;



      void setAdvectionSpeedY(const double advectionSpeed)

      void setAdvectionSpeedY(const Real& advectionSpeed)

      {

	   this->advectionSpeedY = advectionSpeed;

      }





      void setAdvectionSpeedX(const double advectionSpeed)

      void setAdvectionSpeedX(const Real& advectionSpeed)

      {

	   this->advectionSpeedX = advectionSpeed;

      }



      void setViscosity(const double artificalViscosity)

      void setViscosity(const Real& artificalViscosity)

      {

	   this->artificalViscosity = artificalViscosity;

      }
@@ -161,22 +161,22 @@ class LaxFridrichs< tnlGrid< 3,MeshReal, Device, MeshIndex >, Real, Index > 
      typedef tnlMeshFunction< MeshType > MeshFunctionType;

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

      Real tau;

      double artificalViscosity;

      double advectionSpeedX;

      double advectionSpeedY;

      Real artificalViscosity;

      Real advectionSpeedX;

      Real advectionSpeedY;



      void setAdvectionSpeedY(const double advectionSpeed)

      void setAdvectionSpeedY(const Real& advectionSpeed)

      {

	   this->advectionSpeedY = advectionSpeed;

      }





      void setAdvectionSpeedX(const double advectionSpeed)

      void setAdvectionSpeedX(const Real& advectionSpeed)

      {

	   this->advectionSpeedX = advectionSpeed;

      }



      void setViscosity(const double artificalViscosity)

      void setViscosity(const Real& artificalViscosity)

      {

	   this->artificalViscosity = artificalViscosity;

      }

examples/advection/LaxFridrichs_impl.h

+6 −18
Original line number Diff line number Diff line
@@ -48,12 +48,8 @@ operator()( const MeshFunction& u, 
   double a;

   a = this->advectionSpeedX;

   return   (0.5 / this->tau ) * this->artificalViscosity *

	    ( u[ west ]

            - 2.0 * u[ center ]

            + u[ east ] )

            - (a * ( u[ east ]

            - u[west] ) )

            * hxInverse * 0.5;

	    ( u[ west ]- 2.0 * u[ center ] + u[ east ] )

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



}
@@ -165,18 +161,10 @@ operator()( const MeshFunction& u, 
   double b;

   a = this->advectionSpeedX;

   b = this->advectionSpeedY;

   return ( 0.25 / this->tau ) * this->artificalViscosity * (

              u[ west ]

            + u[ east ]

            + u[ south ]

            + u[ north ]

            - 4 * u[ center ] )

            - (a * ( u[ east ]

            - u[west] ) )

            * hxInverse * 0.5

            - (b * ( u[ south ]

            - u[ south ] ) )

            * hyInverse * 0.5;

   return ( 0.25 / this->tau ) * this->artificalViscosity * 

          ( u[ west ] + u[ east ] + u[ south ] + u[ north ] - 4 * u[ center ] ) -

          (a * ( u[ east ] - u[west] ) ) * hxInverse * 0.5 - 

          (b * ( u[ north ] - u[ south ] ) ) * hyInverse * 0.5;





}

examples/advection/advection.h

+4 −0
Original line number Diff line number Diff line
@@ -40,6 +40,10 @@ template< typename ConfigTag >class advectionConfig 
	 config.addEntry< tnlString >( "move", "choose movement type", "advection");

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

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

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

	    config.addEntryEnum< int >( 1 );

	    config.addEntryEnum< int >( 2 );

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



         /****

          * Add definition of your solver command line arguments.

examples/advection/advectionProblem_impl.h

+39 −31
Original line number Diff line number Diff line
@@ -94,35 +94,40 @@ setInitialCondition( const tnlParameterContainer& parameters, 
                     DofVectorType& dofs,

                     MeshDependentDataType& meshDependentData )

{

   cout << "vaules adding";

   typedef typename MeshType::Cell Cell;

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

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

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

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

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

   int dimensions = 2; //vyresit!!!!

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

   getchar();

   if (beginChoice == "sin_square")

      {

	   double constantFunction;

	   if (dimensions == 1)

	       {

                   cout << "adding DOFS" << endl;

		   dofs[0] = 0;

		   double expValue;

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

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

		   {

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

			if ((i>0.4*count) && (i<0.5*count)) constantFunction=1; else constantFunction=0;

			if (expValue>constantFunction) dofs[i] = expValue; else dofs[i] = constantFunction;

		   };

		   dofs[count] = 0;

		};

	    if (dimensions == 2)

		   dofs[count-1] = 0;

		}

	    else if (dimensions == 2)

	       {

                   count = sqrt(count);

		   double expValue;

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

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

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

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

		      {

			expValue = exp(-pow(size*i-2,2)-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-1) * count -1+ j] = expValue; else dofs[(i-1) * count -1+ j] = constantFunction;

			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;

		      };

		};

       }
@@ -131,24 +136,25 @@ setInitialCondition( const tnlParameterContainer& parameters, 
	   if (dimensions == 1)

	      {

		   dofs[0] = 0;

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

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

		   {

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

		   };

		   dofs[count] = 0;

		};

	    if (dimensions == 2)

		   dofs[count-1] = 0;

		}

	    else if (dimensions == 2)

	       {

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

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

                   count = sqrt(count);

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

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

		      {

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

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

		      };

		};

     };

   //setting velocity field



   

   cout << dofs << endl;

   getchar();

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

   if( ! dofs.load( initialConditionFile ) )

   {
@@ -164,6 +170,7 @@ setInitialCondition( const tnlParameterContainer& parameters, 
   differentialOperator.setAdvectionSpeedX(advectionSpeedX);

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

   differentialOperator.setAdvectionSpeedY(advectionSpeedY);

   cout << "vaules added";

   return true;

}
@@ -236,9 +243,9 @@ getExplicitRHS( const RealType& time, 
    * You may use supporting mesh dependent data if you need.

    */

   typedef typename MeshType::Cell Cell;

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

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

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

   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++)
@@ -259,11 +266,12 @@ getExplicitRHS( const RealType& time, 
            ;}

  }

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

  { 

*/  { 

   this->bindDofs( mesh, _u );

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

   MeshFunctionType u( mesh, _u ); 

   MeshFunctionType fu( mesh, _fu );

   differentialOperator.setTau(tau); 

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

                                                           mesh,

                                                           this->differentialOperator,
@@ -271,11 +279,11 @@ getExplicitRHS( const RealType& time, 
                                                           this->rightHandSide,

                                                           u,

                                                           fu );

   tnlBoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter; 

/*   tnlBoundaryConditionsSetter< MeshFunctionType, BoundaryCondition > boundaryConditionsSetter; 

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

      this->boundaryCondition, 

      time + tau, 

       u ); 

       u ); */

 }

}

template< typename Mesh,

examples/advection/tnl-run-advection

+2 −2
Original line number Diff line number Diff line
@@ -7,10 +7,10 @@ tnl-grid-setup --dimensions 1 \ 
 

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

#         --output-file init.tnl

./advection --time-discretisation explicit \

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

	      --time-step 1.0e-3 \

              --boundary-conditions-constant 0 \

              --discrete-solver euler \

              --discrete-solver merson \

              --snapshot-period 0.1 \

              --final-time 1.0 \

	      --artifical-viscosity 1.0 \