Skip to content
Snippets Groups Projects
Commit 7affe081 authored by Jan Schäfer's avatar Jan Schäfer
Browse files

dodelana analyticka cast advekce

rozctvrceni pocatecni podminky u eulera 2D
parent 153def2c
No related branches found
No related tags found
No related merge requests found
...@@ -341,7 +341,7 @@ getExplicitRHS( const RealType& time, ...@@ -341,7 +341,7 @@ getExplicitRHS( const RealType& time,
RealType expValue; RealType expValue;
for (IndexType i = 1; i < count-2; i++) for (IndexType i = 1; i < count-2; i++)
{ {
expValue = exp(-pow(this->size*i-0.2*size-step * tau * this->speedX,2)); expValue = exp(-pow(10/(size*count)*(size*i-0.2*size*count)-step * 10 * 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 ((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; if (expValue>constantFunction) this->analyt[i] = expValue; else this->analyt[i] = constantFunction;
}; };
...@@ -353,7 +353,7 @@ getExplicitRHS( const RealType& time, ...@@ -353,7 +353,7 @@ getExplicitRHS( const RealType& time,
for (IndexType i = 0; i < inverseSquareCount-1; i++) for (IndexType i = 0; i < inverseSquareCount-1; i++)
for (IndexType j = 0; j < inverseSquareCount-1; j++) 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)); expValue = 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));
if ((i - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedX>0.4*inverseSquareCount) && (i - step * tau * (inverseSquareCount/this->schemeSize) * this -> speedX<0.5*inverseSquareCount) 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)) && (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; constantFunction=1; else constantFunction=0;
...@@ -369,7 +369,7 @@ getExplicitRHS( const RealType& time, ...@@ -369,7 +369,7 @@ getExplicitRHS( const RealType& time,
this->analyt[0] = 0; this->analyt[0] = 0;
for (IndexType i = 1; i < count-2; i++) 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[i] = exp(-pow(10/(size*count)*(size*i-0.2*size*count)-step * 10 * tau * this->speedX,2));
}; };
this->analyt[count-1] = 0; this->analyt[count-1] = 0;
} }
...@@ -379,7 +379,7 @@ getExplicitRHS( const RealType& time, ...@@ -379,7 +379,7 @@ getExplicitRHS( const RealType& time,
for (IndexType i = 1; i < inverseSquareCount-1; i++) for (IndexType i = 1; i < inverseSquareCount-1; i++)
for (IndexType j = 1; j < inverseSquareCount-1; j++) 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)); 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));
}; };
}; };
} }
......
...@@ -33,15 +33,23 @@ template< typename ConfigTag >class eulerConfig ...@@ -33,15 +33,23 @@ template< typename ConfigTag >class eulerConfig
config.addEntryEnum< tnlString >( "neumann" ); config.addEntryEnum< tnlString >( "neumann" );
config.addEntryEnum< tnlString >( "mymixed" ); config.addEntryEnum< tnlString >( "mymixed" );
config.addEntry< double >( "boundary-conditions-constant", "This sets a value in case of the constant boundary conditions." ); 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-up-density", "This sets a value of left up density." );
config.addEntry< double >( "left-velocityX", "This sets a value of left_x velocity." ); config.addEntry< double >( "left-up-velocityX", "This sets a value of left up x velocity." );
config.addEntry< double >( "left-velocityY", "This sets a value of left_y velocity." ); config.addEntry< double >( "left-up-velocityY", "This sets a value of left up y velocity." );
config.addEntry< double >( "left-pressure", "This sets a value of left pressure." ); config.addEntry< double >( "left-up-pressure", "This sets a value of left up pressure." );
config.addEntry< double >( "riemann-border", "This sets a position of discontinuity." ); config.addEntry< double >( "left-down-density", "This sets a value of left down density." );
config.addEntry< double >( "right-density", "This sets a value of right density." ); config.addEntry< double >( "left-down-velocityX", "This sets a value of left down x velocity." );
config.addEntry< double >( "right-velocityX", "This sets a value of right_x velocity." ); config.addEntry< double >( "left-down-velocityY", "This sets a value of left down y velocity." );
config.addEntry< double >( "right-velocityY", "This sets a value of right_y velocity." ); config.addEntry< double >( "left-down-pressure", "This sets a value of left down pressure." );
config.addEntry< double >( "right-pressure", "This sets a value of right pressure." ); config.addEntry< double >( "riemann-border", "This sets a position of discontinuity cross." );
config.addEntry< double >( "right-up-density", "This sets a value of right density." );
config.addEntry< double >( "right-up-velocityX", "This sets a value of right_x velocity." );
config.addEntry< double >( "right-up-velocityY", "This sets a value of right_y velocity." );
config.addEntry< double >( "right-up-pressure", "This sets a value of right pressure." );
config.addEntry< double >( "right-down-density", "This sets a value of right density." );
config.addEntry< double >( "right-down-velocityX", "This sets a value of right_x velocity." );
config.addEntry< double >( "right-down-velocityY", "This sets a value of right_y velocity." );
config.addEntry< double >( "right-down-pressure", "This sets a value of right pressure." );
config.addEntry< double >( "gamma", "This sets a value of gamma constant." ); config.addEntry< double >( "gamma", "This sets a value of gamma constant." );
/**** /****
......
...@@ -104,16 +104,26 @@ setInitialCondition( const tnlParameterContainer& parameters, ...@@ -104,16 +104,26 @@ setInitialCondition( const tnlParameterContainer& parameters,
{ {
typedef typename MeshType::Cell Cell; typedef typename MeshType::Cell Cell;
gamma = parameters.getParameter< RealType >( "gamma" ); gamma = parameters.getParameter< RealType >( "gamma" );
RealType rhoL = parameters.getParameter< RealType >( "left-density" ); RealType rhoLu = parameters.getParameter< RealType >( "left-up-density" );
RealType velLX = parameters.getParameter< RealType >( "left-velocityX" ); RealType velLuX = parameters.getParameter< RealType >( "left-up-velocityX" );
RealType velLY = parameters.getParameter< RealType >( "left-velocityY" ); RealType velLuY = parameters.getParameter< RealType >( "left-up-velocityY" );
RealType preL = parameters.getParameter< RealType >( "left-pressure" ); RealType preLu = parameters.getParameter< RealType >( "left-up-pressure" );
RealType eL = ( preL / (gamma - 1) ) + 0.5 * rhoL * pow(velLX,2)+pow(velLY,2); RealType eLu = ( preLu / (gamma - 1) ) + 0.5 * rhoLu * pow(velLuX,2)+pow(velLuY,2);
RealType rhoR = parameters.getParameter< RealType >( "right-density" ); RealType rhoLd = parameters.getParameter< RealType >( "left-down-density" );
RealType velRX = parameters.getParameter< RealType >( "right-velocityX" ); RealType velLdX = parameters.getParameter< RealType >( "left-down-velocityX" );
RealType velRY = parameters.getParameter< RealType >( "right-velocityY" ); RealType velLdY = parameters.getParameter< RealType >( "left-down-velocityY" );
RealType preR = parameters.getParameter< RealType >( "right-pressure" ); RealType preLd = parameters.getParameter< RealType >( "left-down-pressure" );
RealType eR = ( preR / (gamma - 1) ) + 0.5 * rhoR * pow(velRX,2)+pow(velRY,2); RealType eLd = ( preLd / (gamma - 1) ) + 0.5 * rhoLd * pow(velLdX,2)+pow(velLdY,2);
RealType rhoRu = parameters.getParameter< RealType >( "right-up-density" );
RealType velRuX = parameters.getParameter< RealType >( "right-up-velocityX" );
RealType velRuY = parameters.getParameter< RealType >( "right-up-velocityY" );
RealType preRu = parameters.getParameter< RealType >( "right-up-pressure" );
RealType eRu = ( preRu / (gamma - 1) ) + 0.5 * rhoRu * pow(velRuX,2)+pow(velRuY,2);
RealType rhoRd = parameters.getParameter< RealType >( "right-down-density" );
RealType velRdX = parameters.getParameter< RealType >( "right-down-velocityX" );
RealType velRdY = parameters.getParameter< RealType >( "right-down-velocityY" );
RealType preRd = parameters.getParameter< RealType >( "right-down-pressure" );
RealType eRd = ( preRd / (gamma - 1) ) + 0.5 * rhoRd * pow(velRdX,2)+pow(velRdY,2);
RealType x0 = parameters.getParameter< RealType >( "riemann-border" ); RealType x0 = parameters.getParameter< RealType >( "riemann-border" );
int size = mesh.template getEntitiesCount< Cell >(); int size = mesh.template getEntitiesCount< Cell >();
uRho.bind(mesh, dofs, 0); uRho.bind(mesh, dofs, 0);
...@@ -130,25 +140,49 @@ setInitialCondition( const tnlParameterContainer& parameters, ...@@ -130,25 +140,49 @@ setInitialCondition( const tnlParameterContainer& parameters,
for(IndexType i = 0; i < sqrt(size); i++) for(IndexType i = 0; i < sqrt(size); i++)
if ((i < x0 * sqrt(size))&&(j < x0 * sqrt(size)) ) if ((i < x0 * sqrt(size))&&(j < x0 * sqrt(size)) )
{ {
uRho[j*sqrt(size)+i] = rhoL; uRho[j*sqrt(size)+i] = rhoLd;
uRhoVelocityX[j*sqrt(size)+i] = rhoL * velLX; uRhoVelocityX[j*sqrt(size)+i] = rhoLd * velLdX;
uRhoVelocityY[j*sqrt(size)+i] = rhoL * velLY; uRhoVelocityY[j*sqrt(size)+i] = rhoLd * velLdY;
uEnergy[j*sqrt(size)+i] = eL; uEnergy[j*sqrt(size)+i] = eL;
velocity[j*sqrt(size)+i] = sqrt(pow(velLX,2)+pow(velLY,2)); velocity[j*sqrt(size)+i] = sqrt(pow(velLdX,2)+pow(velLdY,2));
velocityX[j*sqrt(size)+i] = velLX; velocityX[j*sqrt(size)+i] = velLdX;
velocityY[j*sqrt(size)+i] = velLY; velocityY[j*sqrt(size)+i] = velLdY;
pressure[j*sqrt(size)+i] = preL; pressure[j*sqrt(size)+i] = preLd;
} }
else else
if ((i >= x0 * sqrt(size))&&(j < x0 * sqrt(size)) )
{
uRho[j*sqrt(size)+i] = rhoLu;
uRhoVelocityX[j*sqrt(size)+i] = rhoLu * velLXu;
uRhoVelocityY[j*sqrt(size)+i] = rhoLu * velLYu;
uEnergy[j*sqrt(size)+i] = eLu;
velocity[j*sqrt(size)+i] = sqrt(pow(velLuX,2)+pow(velLuY,2));
velocityX[j*sqrt(size)+i] = velLuX;
velocityY[j*sqrt(size)+i] = velLuY;
pressure[j*sqrt(size)+i] = preLu;
}
else
if ((i >= x0 * sqrt(size))&&(j >= x0 * sqrt(size)) )
{ {
uRho[j*sqrt(size)+i] = rhoR; uRho[j*sqrt(size)+i] = rhoR;
uRhoVelocityX[j*sqrt(size)+i] = rhoR * velRX; uRhoVelocityX[j*sqrt(size)+i] = rhoRu * velRuX;
uRhoVelocityY[j*sqrt(size)+i] = rhoR * velRY; uRhoVelocityY[j*sqrt(size)+i] = rhoRu * velRuY;
uEnergy[j*sqrt(size)+i] = eR; uEnergy[j*sqrt(size)+i] = eRu;
velocity[j*sqrt(size)+i] = sqrt(pow(velRX,2)+pow(velRY,2)); velocity[j*sqrt(size)+i] = sqrt(pow(velRuX,2)+pow(velRuY,2));
velocityX[j*sqrt(size)+i] = velRX; velocityX[j*sqrt(size)+i] = velRuX;
velocityY[j*sqrt(size)+i] = velRY; velocityY[j*sqrt(size)+i] = velRuY;
pressure[j*sqrt(size)+i] = preR; pressure[j*sqrt(size)+i] = preRu;
}
else
{
uRho[j*sqrt(size)+i] = rhoRd;
uRhoVelocityX[j*sqrt(size)+i] = rhoRd * velRdX;
uRhoVelocityY[j*sqrt(size)+i] = rhoRd * velRdY;
uEnergy[j*sqrt(size)+i] = eRd;
velocity[j*sqrt(size)+i] = sqrt(pow(velRdX,2)+pow(velRdY,2));
velocityX[j*sqrt(size)+i] = velRdX;
velocityY[j*sqrt(size)+i] = velRdY;
pressure[j*sqrt(size)+i] = preRd;
}; };
return true; return true;
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment