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

zobecneny pocatecni podminky

parent 046267bd

examples/advection/advection.h

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

#include "advectionRhs.h"

#include "advectionBuildConfigTag.h"

#include "tnlRiemann1DBoundaryConditions.h"

#include "tnlRiemann2DBoundaryConditions.h"



typedef advectionBuildConfigTag BuildConfig;
@@ -30,11 +32,15 @@ template< typename ConfigTag >class advectionConfig 
         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 >( "sin");

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

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

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

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

	    config.addEntryEnum< tnlString >( "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< tnlString >( "move", "choose movement type", "advection");
@@ -102,11 +108,28 @@ class advectionSetter 
             SolverStarter solverStarter;

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

          }

          if( boundaryConditionsType == "riemann1D" )

          {

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

          }

      }



};

examples/advection/advectionProblem.h

+8 −0
Original line number Diff line number Diff line
@@ -26,6 +26,14 @@ class advectionProblem: 
      using typename BaseType::DofVectorType;

      using typename BaseType::MeshDependentDataType;

      static tnlString getTypeStatic();

      int dimension;

      tnlString choice;

      RealType size;

      int step = 0;

      MeshFunctionType analyt;

      RealType speedX;

      RealType speedY;

      RealType schemeSize;



      tnlString getPrologHeader() const;

examples/advection/advectionProblem_impl.h

+176 −63
Original line number Diff line number Diff line
@@ -94,26 +94,50 @@ 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 >();

   int inveseSquareCount = sqrt(count);

   int inverseSquareCount = sqrt(count);

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

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

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

   getchar();

   if (beginChoice == "sin_square")

   tnlVector < RealType, DeviceType, IndexType > data1;

   data1.setSize(count);

   this->analyt.bind(mesh,data1);

   this->dimension = dimensions;

   this->choice = beginChoice;

   this->size = size;

   this->schemeSize =parameters.getParameter< RealType >( "realSize" );

   if (beginChoice == "square")

      {

	   if (dimensions == 1)

	       {

		   dofs[0] = 0;

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

		   {

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

		   };

		   dofs[count-1] = 0;

		}

	    else if (dimensions == 2)

	       {

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

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

		      {

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

                        dofs[i * inverseSquareCount + j]=1; else dofs[i * inverseSquareCount + j]=0;

		      };

		};

       }

   else if (beginChoice == "exp_square")

      {

	   RealType constantFunction;

	   if (dimensions == 1)

	       {

                   cout << "adding DOFS" << endl;

		   dofs[0] = 0;

		   RealType expValue;

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

		   {

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

			expValue = exp(-pow(10/(size*count)*(size*i-0.2*size*count),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;

		   };
@@ -122,42 +146,60 @@ setInitialCondition( const tnlParameterContainer& parameters, 
	    else if (dimensions == 2)

	       {

		   RealType expValue;

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

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

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

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

		      {

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

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

			expValue = exp(-pow(10/(size*inverseSquareCount)*(size*i-0.2*size*inverseSquareCount),2)-pow(10/(size*inverseSquareCount)*(size*j-0.2*size*inverseSquareCount),2));

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

                        constantFunction=1; else constantFunction=0;

			if (expValue>constantFunction) dofs[i * inveseSquareCount + j] = expValue; else dofs[i * inveseSquareCount + j]

			if (expValue>constantFunction) dofs[i * inverseSquareCount + j] = expValue; else dofs[i * inverseSquareCount + j]

                         = constantFunction;

		      };

		};

       }

   else if (beginChoice == "sin")

   else if (beginChoice == "exp")

      {

	   if (dimensions == 1)

	      {

		   dofs[0] = 0;

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

		   {

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

			dofs[i] = exp(-pow(10/(size*count)*(size*i-0.2*size*count),2));

		   };

		   dofs[count-1] = 0;

		}

	    else if (dimensions == 2)

	       {

                   count = sqrt(count);

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

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

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

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

		      {

			   dofs[i * inverseSquareCount + j] = exp(-pow(10/(size*inverseSquareCount)*(size*i-0.2*size*inverseSquareCount),2)-pow(10/(size*inverseSquareCount)*(size*j-0.2*size*inverseSquareCount)*size,2));

		      };

		};

     }

   else if (beginChoice == "riemann")

      {

	   if (dimensions == 1)

	      {

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

		   dofs[0] = 0;

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

		   {

			if (i < 0.5 * count ) dofs[i] = 1; else

                        dofs[i] = 0;

		   };

		   dofs[count-1] = 0;

		}

	    else if (dimensions == 2)

	       {

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

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

		      {

			   if (i < 0.5*inverseSquareCount && j < 0.5*inverseSquareCount) dofs[i * inverseSquareCount + j] = 1; else

                           dofs[i * inverseSquareCount + j] = 0;

		      };

		};

     };

   //setting velocity field   

   cout << dofs << endl;

   

   getchar();

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

   if( ! dofs.load( initialConditionFile ) )

   {
@@ -175,6 +217,8 @@ setInitialCondition( const tnlParameterContainer& parameters, 
   velocityY.bind(mesh, data, count);

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

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

   this->speedX =advectionSpeedX;

   this->speedY =advectionSpeedY;

   if (velocityType == "advection")

   {

      for (IndexType i = 0; i<count; i++)
@@ -188,18 +232,14 @@ setInitialCondition( const tnlParameterContainer& parameters, 
      RealType radius;

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

         {

            radius = sqrt(pow(i/inveseSquareCount - (inveseSquareCount/2.0), 2) + pow(i%inveseSquareCount - (inveseSquareCount/2.0), 2))/inveseSquareCount;

            radius = sqrt(pow(i/inverseSquareCount - (inverseSquareCount/2.0), 2) + pow(i%inverseSquareCount - (inverseSquareCount/2.0), 2))/inverseSquareCount;

            if ( radius == 0 ) {velocityX[i] = 0; velocityY[i] = 0; } else

            velocityX[i] = advectionSpeedX * radius * sin(atan(1) * 4 * (i/inveseSquareCount-inveseSquareCount/2.0) / inveseSquareCount);

            velocityY[i] = advectionSpeedX * radius * sin( (-1) * atan(1) * 4 * (i%inveseSquareCount-inveseSquareCount/2.0) / inveseSquareCount);

            velocityX[i] = advectionSpeedX * radius * sin(atan(1) * 4 * (i/inverseSquareCount-inverseSquareCount/2.0) / inverseSquareCount);

            velocityY[i] = advectionSpeedX * radius * sin( (-1) * atan(1) * 4 * (i%inverseSquareCount-inverseSquareCount/2.0) / inverseSquareCount);

         };

   };

   velocityX.write("Vx", "gnuplot");

   velocityY.write("Vy", "gnuplot");

   getchar();

   differentialOperator.setAdvectionSpeedX(velocityX);      

   differentialOperator.setAdvectionSpeedY(velocityY);

   cout << "vaules added";

   return true;

}
@@ -247,6 +287,9 @@ makeSnapshot( const RealType& time, 
   FileNameBaseNumberEnding( "u-", step, 5, ".tnl", fileName );

   if( ! dofs.save( fileName ) )

      return false;

   FileNameBaseNumberEnding( "a-", step, 5, ".tnl", fileName );

   if( ! this->analyt.save( fileName ) )

      return false;

   return true;

}
@@ -263,39 +306,109 @@ getExplicitRHS( const RealType& time, 
                DofVectorType& _fu,

                MeshDependentDataType& meshDependentData )

{

   /****

    * 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.

    */

//   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]))

            ;}

   step++;

   typedef typename MeshType::Cell Cell;

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

   double inverseSquareCount = sqrt(count);

   if (this->choice == "square")

      {

	   if (dimension == 1)

	       {

		   this->analyt[0];

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

		   {

			if ((i - step * tau * (count/this->schemeSize) * this -> speedX>0.4*count) && (i - step * tau * (count/this->schemeSize) * this -> speedX<0.5*count)) analyt[i]=1; else analyt[i]=0;

		   };

		   this->analyt[count-1] = 0;

		}

	    else if (dimension == 2)

	       {

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

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

		      {

			if ((i - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedX>0.4*inverseSquareCount) && (i - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedX<0.5*inverseSquareCount) 

                         && (j - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedY>0.4*inverseSquareCount) && (j - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedY<0.5*inverseSquareCount))

                        analyt[i]=1; else analyt[i]=0;

		      };

		};

       }

   else if (this->choice == "exp_square")

      {

	   RealType constantFunction;

	   if (dimension == 1)

	       {

		   this->analyt[0] = 0;

		   RealType expValue;

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

		   {

			expValue = exp(-pow(this->size*i-0.2*size-step * tau * this->speedX,2));

			if ((i - step * tau * (count/this->schemeSize) * this -> speedX>0.4*count) && (i - step * tau * (count/this->schemeSize) * this -> speedX<0.5*count)) constantFunction=1; else constantFunction=0;

			if (expValue>constantFunction) this->analyt[i] = expValue; else this->analyt[i] = constantFunction;

		   };

		   this->analyt[count-1] = 0;

		}

	    else if (dimension == 2)

	       {

		   RealType expValue;

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

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

		      {

			expValue = exp(-pow(this->size*i-0.2*size-step * tau * this->speedX,2)-pow(this->size*j-0.2*size-step * tau * this->speedY,2));

			if ((i - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedX>0.4*inverseSquareCount) && (i - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedX<0.5*inverseSquareCount) 

                         && (j - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedY>0.4*inverseSquareCount) && (j - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedY<0.5*inverseSquareCount))

                        constantFunction=1; else constantFunction=0;

			if (expValue>constantFunction) this->analyt[i * inverseSquareCount + j] = expValue; else this->analyt[i * inverseSquareCount + j]

                         = constantFunction;

		      };

		};

       }

   else if (this->choice == "exp")

      {

	   if (dimension == 1)

	      {

		   this->analyt[0] = 0;

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

		   {

			this->analyt[i] = exp(-pow(this->size*i-0.2*size-step * tau * this->speedX,2));

		   };

		   this->analyt[count-1] = 0;

		}

	    else if (dimension == 2)

	       {

                   count = sqrt(count);

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

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

		      {

			   this->analyt[i * inverseSquareCount + j] = exp(-pow(this->size*i-0.2*size-step * tau * this->speedX,2)-pow(this->size*j-0.2*size-step * tau * this->speedY,2));

		      };

		};

     }

   else if (this->choice == "riemann")

      {

	   if (dimension == 1)

	      {

		   this->analyt[0] = 0;

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

		   {

			if (i - step * tau * (count/this->schemeSize) * this -> speedX < 0.5 * count ) this->analyt[i] = 1; else

                        this->analyt[i] = 0;

		   };

		   this->analyt[count-1] = 0;

		}

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

*/  { 

	    else if (dimension == 2)

	       {

                   count = sqrt(count);

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

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

		      {

			   if (i - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedX < 0.5*inverseSquareCount && 

                               j - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedY < 0.5*inverseSquareCount) 

                              this->analyt[i * inverseSquareCount + j] = 1; 

                           else

                              this->analyt[i * inverseSquareCount + j] = 0;

		      };

		};

     };

   this->bindDofs( mesh, _u );

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

   MeshFunctionType u( mesh, _u );
@@ -313,7 +426,7 @@ getExplicitRHS( const RealType& time, 
      this->boundaryCondition, 

      time + tau, 

       u ); */

 }



}

template< typename Mesh,

          typename BoundaryCondition,

examples/advection/tnl-run-advection

+1 −1
Original line number Diff line number Diff line
@@ -14,7 +14,7 @@ tnl-advection-dbg --time-discretisation explicit \ 
              --snapshot-period 0.1 \

              --final-time 1.0 \

	      --artifical-viscosity 1.0 \

	      --begin sin \

	      --begin exp \

	      --advection-speedX 2.0 \



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

examples/advection/tnl-run-advection1

+1 −1
Original line number Diff line number Diff line
@@ -17,7 +17,7 @@ tnl-grid-setup --dimensions 2 \ 
              --snapshot-period 0.1 \

              --final-time 1.0 \

	      --artifical-viscosity 0.2 \

	      --begin sin_square \

	      --begin exp_square \

	      --advection-speedX 2.0 \

	      --advection-speedY 2.0 \