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

oprava analyt reseni advekce 

pridani vlastnich neumannovych op
parent afc7f13d

examples/advection/advectionProblem_impl.h

+10 −9
Original line number Diff line number Diff line
@@ -173,7 +173,7 @@ setInitialCondition( const tnlParameterContainer& parameters, 
		   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));

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

		      };

		};

     }
@@ -284,13 +284,14 @@ makeSnapshot( const RealType& time, 
   cout << endl << "Writing output at time " << time << " step " << step << "." << endl;

   this->bindDofs( mesh, dofs );

   tnlString fileName;

   MeshFunctionType dofsh;

   dofsh.bind(mesh,dofs);

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

   if( ! dofs.save( fileName ) )

   if( ! dofsh.save( fileName ) )

      return false;

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

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

      return false;

cin.ignore();

   return true;

}
@@ -317,7 +318,7 @@ getExplicitRHS( const RealType& time, 
      {

	   if (dimension == 1)

	       {

		   this->analyt[0];

		   this->analyt[0] = 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;
@@ -331,7 +332,7 @@ getExplicitRHS( const RealType& time, 
		      {

			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;

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

		      };

		};

       }
@@ -379,8 +380,8 @@ getExplicitRHS( const RealType& time, 
	    else if (dimension == 2)

	       {

                   count = sqrt(count);

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

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

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

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

		      {

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

		      };
@@ -401,8 +402,8 @@ getExplicitRHS( const RealType& time, 
	    else if (dimension == 2)

	       {

                   count = sqrt(count);

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

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

		   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.5*inverseSquareCount && 

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

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

+3 −1
Original line number Diff line number Diff line
@@ -39,7 +39,9 @@ class EulerPressureGetter 
      __cuda_callable__

      Real operator[]( const IndexType& idx ) const

      {

          if (this->rho[ idx ]==0) return 0; else return ( this->gamma - 1.0 ) * ( this->energy[ idx ] - 0.5 * this->rhoVel[ idx ] * this->rhoVel[ idx ] / this->rho[ idx ]);

          //if (this->rho[ idx ]==0) return 0; else return ( this->gamma - 1.0 ) * ( this->energy[ idx ] - 0.5 * this->rhoVel[ idx ] * this->rhoVel[ idx ] / this->rho[ idx ]);

          return ( this->gamma - 1.0 ) * this->energy[ idx ] * this->rho[ idx ];



      }

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

+9 −0
Original line number Diff line number Diff line
@@ -7,6 +7,7 @@ 
#include "eulerProblem.h"

#include "LaxFridrichs.h"

#include "tnlMyMixedBoundaryConditions.h"

#include "tnlMyNeumannBoundaryConditions.h"



#include "eulerRhs.h"

#include "eulerBuildConfigTag.h"
@@ -33,6 +34,7 @@ template< typename ConfigTag >class eulerConfig 
            config.addEntryEnum< tnlString >( "dirichlet" );

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

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

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

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

         config.addEntry< double >( "left-density", "This sets a value of left density." );

         config.addEntry< double >( "left-velocity", "This sets a value of left velocity." );
@@ -107,6 +109,13 @@ class eulerSetter 
             SolverStarter solverStarter;

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

          }

          if( boundaryConditionsType == "myneumann" )

          {

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

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

             SolverStarter solverStarter;

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

          }

          if( boundaryConditionsType == "neumann" )

          {

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

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

+5 −14
Original line number Diff line number Diff line
@@ -99,21 +99,20 @@ setInitialCondition( const tnlParameterContainer& parameters, 
                     DofVectorType& dofs,

                     MeshDependentDataType& meshDependentData )

{

   cout << endl << "get conditions from CML";

   typedef typename MeshType::Cell Cell;

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

   RealType rhoL = parameters.getParameter< RealType >( "left-density" );

   RealType velL = parameters.getParameter< RealType >( "left-velocity" );

   RealType preL = parameters.getParameter< RealType >( "left-pressure" );

   RealType eL = ( preL / (gamma - 1) ) + 0.5 * rhoL * velL * velL;

   RealType eL = ( preL / ( rhoL * (gamma - 1) ) );

   //RealType eL = ( preL / (gamma - 1) ) + 0.5 * rhoL * velL * velL;

   RealType rhoR = parameters.getParameter< RealType >( "right-density" );

   RealType velR = parameters.getParameter< RealType >( "right-velocity" );

   RealType preR = parameters.getParameter< RealType >( "right-pressure" );

   RealType eR = ( preR / (gamma - 1) ) + 0.5 * rhoR * velR * velR;

   RealType eR = ( preR / ( rhoR * (gamma - 1) ) );

   //RealType eR = ( preR / (gamma - 1) ) + 0.5 * rhoR * velR * velR;

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

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

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

cout << count << endl;

   uRho.bind(mesh, dofs, 0);

   uRhoVelocity.bind(mesh, dofs, count);

   uEnergy.bind(mesh, dofs, 2 * count);
@@ -121,7 +120,6 @@ cout << count << endl; 
   data.setSize(2*count);

   velocity.bind( mesh, data, 0);

   pressure.bind( mesh, data, count );  

   cout << endl << "set conditions from CML"<< endl;   

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

      if (i < x0 * count )

         {
@@ -139,13 +137,6 @@ cout << count << endl; 
            velocity[i] = velR;

            pressure[i] = preR;

         };

   cout << "dofs = " << dofs << endl;

   cout << "velocity = " << velocity.getData() << endl;

   cout << "pressure = " << pressure.getData() << endl;



   getchar();

  

   

   /*

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

   if( ! dofs.load( initialConditionFile ) )

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

0 → 100644
+140 −0
Original line number Diff line number Diff line 
/*** coppied and changed

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

                          tnlMyNeumannBoundaryConditions.h  -  description

                             -------------------

    begin                : Nov 17, 2014

    copyright            : (C) 2014 by oberhuber

    email                : tomas.oberhuber@fjfi.cvut.cz

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



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

 *                                                                         *

 *   This program is free software; you can redistribute it and/or modify  *

 *   it under the terms of the GNU General Public License as published by  *

 *   the Free Software Foundation; either version 2 of the License, or     *

 *   (at your option) any later version.                                   *

 *                                                                         *

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



#ifndef TNLMyNeumannBoundaryConditions_H_

#define TNLMyNeumannBoundaryConditions_H_



#include <operators/tnlOperator.h>

#include <functions/tnlConstantFunction.h>

#include <functions/tnlFunctionAdapter.h>



template< typename Mesh,

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

          int MeshEntitiesDimensions = Mesh::getMeshDimensions(),

          typename Real = typename Mesh::RealType,

          typename Index = typename Mesh::IndexType >

class tnlMyNeumannBoundaryConditions

: public tnlOperator< Mesh,

                      MeshBoundaryDomain,

                      MeshEntitiesDimensions,

                      MeshEntitiesDimensions,

                      Real,

                      Index >

{

   public:



      typedef Mesh MeshType;

      typedef Function FunctionType;

      typedef Real RealType;

      typedef typename MeshType::DeviceType DeviceType;

      typedef Index IndexType;



      typedef tnlVector< RealType, DeviceType, IndexType> DofVectorType;

      typedef typename MeshType::VertexType VertexType;



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



      static void configSetup( tnlConfigDescription& config,

                               const tnlString& prefix = "" )

      {

         Function::configSetup( config, prefix );

      }

      

      bool setup( const tnlParameterContainer& parameters,

                  const tnlString& prefix = "" )

      {

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

      }



      void setFunction( const Function& function )

      {

         this->function = function;

      }



      Function& getFunction()

      {

         return this->function;

      }

      

      const Function& getFunction() const

      {

         return this->function;

      }



      template< typename EntityType,

                typename MeshFunction >

      __cuda_callable__

      const RealType operator()( const MeshFunction& u,

                                 const EntityType& entity,                            

                                 const RealType& time = 0 ) const

      {

      const MeshType& mesh = entity.getMesh();

      const auto& neighbourEntities = entity.getNeighbourEntities();

      const IndexType& index = entity.getIndex();

      if( entity.getCoordinates().x() == 0 )

         return u[ neighbourEntities.template getEntityIndex< 1 >() ];

      else

         return u[ neighbourEntities.template getEntityIndex< -1 >() ];  

         //tady se asi delaji okrajove podminky

         //static_assert( EntityType::getDimensions() == MeshEntitiesDimensions, "Wrong mesh entity dimensions." );

      }



      template< typename EntityType >

      __cuda_callable__

      IndexType getLinearSystemRowLength( const MeshType& mesh,

                                          const IndexType& index,

                                          const EntityType& entity ) const

      {

         return 1;

      }



      template< typename PreimageFunction,

                typename MeshEntity,

                typename Matrix,

                typename Vector >

      __cuda_callable__

      void setMatrixElements( const PreimageFunction& u,

                              const MeshEntity& entity,

                              const RealType& time,

                              const RealType& tau,

                              Matrix& matrix,

                              Vector& b ) const

      {

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

      }

   



   protected:



      Function function;

   

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

};



template< typename Mesh,

          typename Function >

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

{

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

   return str;

}



#endif /* TNLMyNeumannBoundaryConditions_H_ */