From 50e651dc73d7b263c186a363b0a55ecc4a81c354 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jan=20Sch=C3=A4fer?= <schafjan@fjfi.cvut.cz>
Date: Wed, 7 Dec 2016 18:04:54 +0100
Subject: [PATCH] opraveny chyby, 1D a 2D ozkouseno
---
.../inviscid-flow/1d/Euler1DPressureGetter.h | 4 +-
.../1d/LaxFridrichsContinuity_impl.h | 2 +
.../1d/LaxFridrichsMomentum_impl.h | 4 +
examples/inviscid-flow/1d/eulerProblem_impl.h | 20 ++---
.../inviscid-flow/2d/EulerPressureGetter.h | 6 +-
.../2d/LaxFridrichsContinuity_impl .h | 2 +-
.../2d/MyNeumannBoundaryConditions.h | 2 +-
examples/inviscid-flow/2d/euler.h | 1 +
examples/inviscid-flow/2d/eulerProblem_impl.h | 28 ++++---
examples/inviscid-flow/3d/euler.h | 8 ++
examples/inviscid-flow/3d/eulerProblem_impl.h | 76 ++++++++++++-------
11 files changed, 99 insertions(+), 54 deletions(-)
diff --git a/examples/inviscid-flow/1d/Euler1DPressureGetter.h b/examples/inviscid-flow/1d/Euler1DPressureGetter.h
index 7900c4b085..57f3ff859d 100644
--- a/examples/inviscid-flow/1d/Euler1DPressureGetter.h
+++ b/examples/inviscid-flow/1d/Euler1DPressureGetter.h
@@ -42,8 +42,8 @@ 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 ]);
- return ( this->gamma - 1.0 ) * this->energy[ 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 ];
}
diff --git a/examples/inviscid-flow/1d/LaxFridrichsContinuity_impl.h b/examples/inviscid-flow/1d/LaxFridrichsContinuity_impl.h
index 3880bb7d2d..dd13d36428 100644
--- a/examples/inviscid-flow/1d/LaxFridrichsContinuity_impl.h
+++ b/examples/inviscid-flow/1d/LaxFridrichsContinuity_impl.h
@@ -45,6 +45,8 @@ operator()( const MeshFunction& u,
const IndexType& west = neighbourEntities.template getEntityIndex< -1 >();
return (0.5 / this->tau) * ( u[ west ] - 2.0 * u[ center ] + u[ east ] )
- 0.5 * hxInverse * ( u[ east ] * this->velocity[ east ] - u[ west ] * this->velocity[ west ] );
+ /*(0.5) * ( u[ west ] + u[ east ] )
+ - 0.5 * hxInverse * this->tau * ( u[ east ] * this->velocity[ east ] - u[ west ] * this->velocity[ west ] );*/
}
template< typename MeshReal,
diff --git a/examples/inviscid-flow/1d/LaxFridrichsMomentum_impl.h b/examples/inviscid-flow/1d/LaxFridrichsMomentum_impl.h
index 217afc4b5f..ebfe407324 100644
--- a/examples/inviscid-flow/1d/LaxFridrichsMomentum_impl.h
+++ b/examples/inviscid-flow/1d/LaxFridrichsMomentum_impl.h
@@ -47,6 +47,10 @@ operator()( const MeshFunction& u,
- 0.5 * hxInverse *
(( u[ east ] * this -> velocity[ east ] + this -> pressure [ east ] )
-( u[ west ] * this -> velocity[ west ] + this -> pressure [ west ] ));
+ /*(0.5) * ( u[ west ] + u[ east ] )
+ - 0.5 * hxInverse * this->tau *
+ (( u[ east ] * this -> velocity[ east ] + this -> pressure [ east ] )
+ -( u[ west ] * this -> velocity[ west ] + this -> pressure [ west ] ));*/
}
template< typename MeshReal,
diff --git a/examples/inviscid-flow/1d/eulerProblem_impl.h b/examples/inviscid-flow/1d/eulerProblem_impl.h
index 1471a46c2d..82e89903d4 100644
--- a/examples/inviscid-flow/1d/eulerProblem_impl.h
+++ b/examples/inviscid-flow/1d/eulerProblem_impl.h
@@ -60,7 +60,7 @@ setup( const MeshPointer& meshPointer,
const Config::ParameterContainer& parameters,
const String& prefix )
{
- if( ! this->boundaryConditionsPointer->setup( meshPointer, parameters, prefix + "boundary-conditions-" ) ||
+ if( //! this->boundaryConditionsPointer->setup( meshPointer, parameters, prefix + "boundary-conditions-" ) ||
! this->rightHandSidePointer->setup( parameters, prefix + "right-hand-side-" ) )
return false;
return true;
@@ -103,19 +103,19 @@ setInitialCondition( const Config::ParameterContainer& parameters,
DofVectorPointer& dofs,
MeshDependentDataPointer& meshDependentData )
{
- std::cout << std::endl << "get conditions from CML";
+ std::cout << std::endl << "get conditions from CMD";
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 / ( rhoL * (gamma - 1) ) );
- //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 / ( rhoR * (gamma - 1) ) );
- //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" );
int count = mesh->template getEntitiesCount< Cell >();
uRho->bind( mesh, *dofs, 0);
@@ -125,7 +125,7 @@ setInitialCondition( const Config::ParameterContainer& parameters,
data.setSize(2*count);
velocity->bind( mesh, data, 0);
pressure->bind( mesh, data, count );
- std::cout << std::endl << "set conditions from CML"<< std::endl;
+ std::cout << std::endl << "set conditions from CMD"<< std::endl;
for(IndexType i = 0; i < count; i++)
if (i < x0 * count )
{
@@ -230,9 +230,9 @@ getExplicitRHS( const RealType& time,
typedef typename MeshType::Cell Cell;
int count = mesh->template getEntitiesCount< Cell >();
//bind _fu
- this->fuRho->bind(mesh, _u, 0);
- this->fuRhoVelocity->bind(mesh, _u, count);
- this->fuEnergy->bind(mesh, _u, 2 * count);
+ this->fuRho->bind(mesh, _fu, 0);
+ this->fuRhoVelocity->bind(mesh, _fu, count);
+ this->fuEnergy->bind(mesh, _fu, 2 * count);
//generating Differential operator object
SharedPointer< Continuity > lF1DContinuity;
diff --git a/examples/inviscid-flow/2d/EulerPressureGetter.h b/examples/inviscid-flow/2d/EulerPressureGetter.h
index 76792cdea9..a251372a77 100644
--- a/examples/inviscid-flow/2d/EulerPressureGetter.h
+++ b/examples/inviscid-flow/2d/EulerPressureGetter.h
@@ -43,10 +43,10 @@ class EulerPressureGetter
__cuda_callable__
Real operator[]( const IndexType& idx ) const
{
-/* if (this->rho[ idx ]==0) return 0;
+ if (this->rho[ idx ]==0) {return 0;}
else return ( this->gamma - 1.0 ) * ( this->energy[ idx ] - 0.5 * this->rho[ idx ] *
- ( this->rhoVelX[ idx ] / this->rho[ idx ] + this->rhoVelY[ idx ] / this->rho[ idx ]) );
-*/ return ( this->gamma - 1.0 ) * ( this->energy[ idx ] * this->rho[ idx ] );
+ ( std::pow(this->rhoVelX[ idx ] / this->rho[ idx ],2) + std::pow(this->rhoVelY[ idx ] / this->rho[ idx ],2) ) );
+//*/ return ( this->gamma - 1.0 ) * ( this->energy[ idx ] * this->rho[ idx ] );
}
diff --git a/examples/inviscid-flow/2d/LaxFridrichsContinuity_impl .h b/examples/inviscid-flow/2d/LaxFridrichsContinuity_impl .h
index 11a7e3adc5..539e574c96 100644
--- a/examples/inviscid-flow/2d/LaxFridrichsContinuity_impl .h
+++ b/examples/inviscid-flow/2d/LaxFridrichsContinuity_impl .h
@@ -46,7 +46,7 @@ operator()( const MeshFunction& u,
const IndexType& center = entity.getIndex();
const IndexType& east = neighbourEntities.template getEntityIndex< 1 >();
const IndexType& west = neighbourEntities.template getEntityIndex< -1 >();
- return ( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse;
+ return 0;/*( u[ west ] - 2.0 * u[ center ] + u[ east ] ) * hxSquareInverse;*/
}
template< typename MeshReal,
diff --git a/examples/inviscid-flow/2d/MyNeumannBoundaryConditions.h b/examples/inviscid-flow/2d/MyNeumannBoundaryConditions.h
index 7987b93160..72285a4ae0 100644
--- a/examples/inviscid-flow/2d/MyNeumannBoundaryConditions.h
+++ b/examples/inviscid-flow/2d/MyNeumannBoundaryConditions.h
@@ -66,7 +66,7 @@ class MyNeumannBoundaryConditions
const Config::ParameterContainer& parameters,
const String& prefix = "" )
{
- return Functions::FunctionAdapter< MeshType, FunctionType >::template setup< MeshPointer >( this->function, meshPointer, parameters, prefix );
+ return true; //Functions::FunctionAdapter< MeshType, FunctionType >::template setup< MeshPointer >( this->function, meshPointer, parameters, prefix );
}
void setFunction( const Function& function )
diff --git a/examples/inviscid-flow/2d/euler.h b/examples/inviscid-flow/2d/euler.h
index acdf88f23e..ef1da7735c 100644
--- a/examples/inviscid-flow/2d/euler.h
+++ b/examples/inviscid-flow/2d/euler.h
@@ -79,6 +79,7 @@ class eulerSetter
static bool run( const Config::ParameterContainer & parameters )
{
+std::cout <<"run";
enum { Dimensions = MeshType::getMeshDimensions() };
typedef LaxFridrichs2D< MeshType, Real, Index > ApproximateOperator;
typedef eulerRhs< MeshType, Real > RightHandSide;
diff --git a/examples/inviscid-flow/2d/eulerProblem_impl.h b/examples/inviscid-flow/2d/eulerProblem_impl.h
index 29b0dd0a34..7c07eafeb0 100644
--- a/examples/inviscid-flow/2d/eulerProblem_impl.h
+++ b/examples/inviscid-flow/2d/eulerProblem_impl.h
@@ -10,6 +10,8 @@
#include "EulerPressureGetter.h"
#include "Euler2DVelXGetter.h"
#include "EulerVelGetter.h"
+#include "MyMixedBoundaryConditions.h"
+#include "MyNeumannBoundaryConditions.h"
namespace TNL {
@@ -59,7 +61,7 @@ setup( const MeshPointer& meshPointer,
const Config::ParameterContainer& parameters,
const String& prefix )
{
- if( ! this->boundaryConditionPointer->setup( meshPointer, parameters, prefix + "boundary-conditions-" ) ||
+ if( //! this->boundaryConditionPointer->setup( meshPointer, parameters, prefix + "boundary-conditions-" ) ||
! this->rightHandSidePointer->setup( parameters, prefix + "right-hand-side-" ) )
return false;
return true;
@@ -76,7 +78,7 @@ getDofs( const MeshPointer& mesh ) const
/****
* Return number of DOFs (degrees of freedom) i.e. number
* of unknowns to be resolved by the main solver.
- */
+ */
return 4*mesh->template getEntitiesCount< typename MeshType::Cell >();
}
@@ -108,26 +110,30 @@ setInitialCondition( const Config::ParameterContainer& parameters,
RealType velLuX = parameters.getParameter< RealType >( "NW-velocityX" );
RealType velLuY = parameters.getParameter< RealType >( "NW-velocityY" );
RealType preLu = parameters.getParameter< RealType >( "NW-pressure" );
- RealType eLu = ( preLu / ( rhoLu * (gamma - 1) ) );
- //RealType eLu = ( preLu / (gamma - 1) ) + 0.5 * rhoLu * pow(velLuX,2)+pow(velLuY,2);
+ //RealType eLu = ( preLu / ( rhoLu * (gamma - 1) ) );
+ RealType eLu = ( preLu / (gamma - 1) ) + 0.5 * rhoLu * ( std::pow(velLuX,2)+std::pow(velLuY,2) );
+ std::cout << rhoLu <<' '<< velLuX<<' '<<velLuY<<' '<<preLu<<' '<<eLu<<' ';
RealType rhoLd = parameters.getParameter< RealType >( "SW-density" );
RealType velLdX = parameters.getParameter< RealType >( "SW-velocityX" );
RealType velLdY = parameters.getParameter< RealType >( "SW-velocityY" );
RealType preLd = parameters.getParameter< RealType >( "SW-pressure" );
- RealType eLd = ( preLd / ( rhoLd * (gamma - 1) ) );
- //RealType eLd = ( preLd / (gamma - 1) ) + 0.5 * rhoLd * pow(velLdX,2)+pow(velLdY,2);
+ //RealType eLd = ( preLd / ( rhoLd * (gamma - 1) ) );
+ RealType eLd = ( preLd / (gamma - 1) ) + 0.5 * rhoLd * ( std::pow(velLdX,2)+std::pow(velLdY,2) );
+ std::cout << rhoLd <<' '<< velLdX<<' '<<velLdY<<' '<<preLd<<' '<<eLd<<' ';
RealType rhoRu = parameters.getParameter< RealType >( "NE-density" );
RealType velRuX = parameters.getParameter< RealType >( "NE-velocityX" );
RealType velRuY = parameters.getParameter< RealType >( "NE-velocityY" );
RealType preRu = parameters.getParameter< RealType >( "NE-pressure" );
- RealType eRu = ( preRu / ( rhoRu * (gamma - 1) ) );
- //RealType eRu = ( preRu / (gamma - 1) ) + 0.5 * rhoRu * pow(velRuX,2)+pow(velRuY,2);
+ //RealType eRu = ( preRu / ( rhoRu * (gamma - 1) ) );
+ RealType eRu = ( preRu / (gamma - 1) ) + 0.5 * rhoRu * ( std::pow(velRuX,2)+std::pow(velRuY,2) );
+ std::cout << rhoRu <<' '<< velRuX<<' '<<velRuY<<' '<<preRu<<' '<<eRu<<' ';
RealType rhoRd = parameters.getParameter< RealType >( "SE-density" );
RealType velRdX = parameters.getParameter< RealType >( "SE-velocityX" );
RealType velRdY = parameters.getParameter< RealType >( "SE-velocityY" );
RealType preRd = parameters.getParameter< RealType >( "SE-pressure" );
- RealType eRd = ( preRd / ( rhoRd * (gamma - 1) ) );
- //RealType eRd = ( preRd / (gamma - 1) ) + 0.5 * rhoRd * pow(velRdX,2)+pow(velRdY,2);
+ //RealType eRd = ( preRd / ( rhoRd * (gamma - 1) ) );
+ RealType eRd = ( preRd / (gamma - 1) ) + 0.5 * rhoRd * ( std::pow(velRdX,2)+std::pow(velRdY,2) );
+ std::cout << rhoRd <<' '<< velRdX<<' '<<velRdY<<' '<<preRd<<' '<<eRd<<' ';
RealType x0 = parameters.getParameter< RealType >( "riemann-border" );
int size = mesh->template getEntitiesCount< Cell >();
uRho->bind(mesh, dofs, 0);
@@ -254,7 +260,7 @@ makeSnapshot( const RealType& time,
fileName.setFileNameBase( "velocity-" );
if( ! velocity->save( fileName.getFileName() ) )
return false;
- fileName.setFileNameBase( "pressue-" );
+ fileName.setFileNameBase( "pressure-" );
if( ! pressure->save( fileName.getFileName() ) )
return false;
diff --git a/examples/inviscid-flow/3d/euler.h b/examples/inviscid-flow/3d/euler.h
index 600653ecc7..d5550d0c72 100644
--- a/examples/inviscid-flow/3d/euler.h
+++ b/examples/inviscid-flow/3d/euler.h
@@ -40,35 +40,43 @@ template< typename ConfigTag >class eulerConfig
config.addEntry< double >( "NWU-density", "This sets a value of northwest up density." );
config.addEntry< double >( "NWU-velocityX", "This sets a value of northwest up x velocity." );
config.addEntry< double >( "NWU-velocityY", "This sets a value of northwest up y velocity." );
+ config.addEntry< double >( "NWU-velocityZ", "This sets a value of northwest up z velocity." );
config.addEntry< double >( "NWU-pressure", "This sets a value of northwest up pressure." );
config.addEntry< double >( "SWU-density", "This sets a value of southwest up density." );
config.addEntry< double >( "SWU-velocityX", "This sets a value of southwest up x velocity." );
config.addEntry< double >( "SWU-velocityY", "This sets a value of southwest up y velocity." );
+ config.addEntry< double >( "SWU-velocityZ", "This sets a value of southwest up z velocity." );
config.addEntry< double >( "SWU-pressure", "This sets a value of southwest up pressure." );
config.addEntry< double >( "NWD-density", "This sets a value of northwest down density." );
config.addEntry< double >( "NWD-velocityX", "This sets a value of northwest down x velocity." );
config.addEntry< double >( "NWD-velocityY", "This sets a value of northwest down y velocity." );
+ config.addEntry< double >( "NWD-velocityZ", "This sets a value of northwest down z velocity." );
config.addEntry< double >( "NWD-pressure", "This sets a value of northwest down pressure." );
config.addEntry< double >( "SWD-density", "This sets a value of southwest down density." );
config.addEntry< double >( "SWD-velocityX", "This sets a value of southwest down x velocity." );
config.addEntry< double >( "SWD-velocityY", "This sets a value of southwest down y velocity." );
+ config.addEntry< double >( "SWF-velocityZ", "This sets a value of southwest down z velocity." );
config.addEntry< double >( "SWD-pressure", "This sets a value of southwest down pressure." );
config.addEntry< double >( "riemann-border", "This sets a position of discontinuity cross." );
config.addEntry< double >( "NEU-density", "This sets a value of northeast up density." );
config.addEntry< double >( "NEU-velocityX", "This sets a value of northeast up x velocity." );
config.addEntry< double >( "NEU-velocityY", "This sets a value of northeast up y velocity." );
+ config.addEntry< double >( "NEU-velocityZ", "This sets a value of northeast up z velocity." );
config.addEntry< double >( "NEU-pressure", "This sets a value of northeast up pressure." );
config.addEntry< double >( "SEU-density", "This sets a value of southeast up density." );
config.addEntry< double >( "SEU-velocityX", "This sets a value of southeast up x velocity." );
config.addEntry< double >( "SEU-velocityY", "This sets a value of southeast up y velocity." );
+ config.addEntry< double >( "SEU-velocityZ", "This sets a value of southeast up z velocity." );
config.addEntry< double >( "SEU-pressure", "This sets a value of southeast up pressure." );
config.addEntry< double >( "NED-density", "This sets a value of northeast down density." );
config.addEntry< double >( "NED-velocityX", "This sets a value of northeast down x velocity." );
config.addEntry< double >( "NED-velocityY", "This sets a value of northeast down y velocity." );
+ config.addEntry< double >( "NED-velocityZ", "This sets a value of northeast down z velocity." );
config.addEntry< double >( "NED-pressure", "This sets a value of northeast down pressure." );
config.addEntry< double >( "SED-density", "This sets a value of southeast down density." );
config.addEntry< double >( "SED-velocityX", "This sets a value of southeast down x velocity." );
config.addEntry< double >( "SED-velocityY", "This sets a value of southeast down y velocity." );
+ config.addEntry< double >( "SED-velocityZ", "This sets a value of southeast down z velocity." );
config.addEntry< double >( "SED-pressure", "This sets a value of southeast down pressure." );
config.addEntry< double >( "gamma", "This sets a value of gamma constant." );
diff --git a/examples/inviscid-flow/3d/eulerProblem_impl.h b/examples/inviscid-flow/3d/eulerProblem_impl.h
index 8f314ec73b..b01979b87c 100644
--- a/examples/inviscid-flow/3d/eulerProblem_impl.h
+++ b/examples/inviscid-flow/3d/eulerProblem_impl.h
@@ -60,7 +60,7 @@ setup( const MeshPointer& meshPointer,
const Config::ParameterContainer& parameters,
const String& prefix )
{
- if( ! this->boundaryConditionPointer->setup( meshPointer, parameters, prefix + "boundary-conditions-" ) ||
+ if( //! this->boundaryConditionPointer->setup( meshPointer, parameters, prefix + "boundary-conditions-" ) ||
! this->rightHandSidePointer->setup( parameters, prefix + "right-hand-side-" ) )
return false;
return true;
@@ -108,51 +108,59 @@ setInitialCondition( const Config::ParameterContainer& parameters,
RealType rhoNWU = parameters.getParameter< RealType >( "NWU-density" );
RealType velNWUX = parameters.getParameter< RealType >( "NWU-velocityX" );
RealType velNWUY = parameters.getParameter< RealType >( "NWU-velocityY" );
+ RealType velNWUZ = parameters.getParameter< RealType >( "NWU-velocityZ" );
RealType preNWU = parameters.getParameter< RealType >( "NWU-pressure" );
- RealType eNWU = ( preNWU / ( rhoNWU * (gamma - 1) ) );
- //RealType eNWU = ( preNWU / (gamma - 1) ) + 0.5 * rhoNWU * pow(velNWUX,2)+pow(velNWUY,2);
+ //RealType eNWU = ( preNWU / ( rhoNWU * (gamma - 1) ) );
+ RealType eNWU = ( preNWU / (gamma - 1) ) + 0.5 * rhoNWU * (std::pow(velNWUX,2)+std::pow(velNWUY,2)+std::pow(velNWUZ,2));
RealType rhoSWU = parameters.getParameter< RealType >( "SWU-density" );
RealType velSWUX = parameters.getParameter< RealType >( "SWU-velocityX" );
RealType velSWUY = parameters.getParameter< RealType >( "SWU-velocityY" );
+ RealType velSWUZ = parameters.getParameter< RealType >( "SWU-velocityZ" );
RealType preSWU = parameters.getParameter< RealType >( "SWU-pressure" );
- RealType eSWU = ( preSWU / ( rhoSWU * (gamma - 1) ) );
- //RealType eSWU = ( preSWU / (gamma - 1) ) + 0.5 * rhoSWU * pow(velSWUX,2)+pow(velSWUY,2);
+ //RealType eSWU = ( preSWU / ( rhoSWU * (gamma - 1) ) );
+ RealType eSWU = ( preSWU / (gamma - 1) ) + 0.5 * rhoSWU * (std::pow(velSWUX,2)+std::pow(velSWUY,2)+std::pow(velSWUZ,2));
RealType rhoNWD = parameters.getParameter< RealType >( "NWD-density" );
RealType velNWDX = parameters.getParameter< RealType >( "NWD-velocityX" );
RealType velNWDY = parameters.getParameter< RealType >( "NWD-velocityY" );
+ RealType velNWDZ = parameters.getParameter< RealType >( "NWD-velocityZ" );
RealType preNWD = parameters.getParameter< RealType >( "NWD-pressure" );
- RealType eNWD = ( preNWD / ( rhoNWD * (gamma - 1) ) );
- //RealType eNWD = ( preNWD / (gamma - 1) ) + 0.5 * rhoNWD * pow(velNWDX,2)+pow(velNWDY,2);
+ //RealType eNWD = ( preNWD / ( rhoNWD * (gamma - 1) ) );
+ RealType eNWD = ( preNWD / (gamma - 1) ) + 0.5 * rhoNWD * (std::pow(velNWDX,2)+std::pow(velNWDY,2)+std::pow(velNWDZ,2));
RealType rhoSWD = parameters.getParameter< RealType >( "SWD-density" );
RealType velSWDX = parameters.getParameter< RealType >( "SWD-velocityX" );
RealType velSWDY = parameters.getParameter< RealType >( "SWD-velocityY" );
+ RealType velSWDZ = parameters.getParameter< RealType >( "SWD-velocityZ" );
RealType preSWD = parameters.getParameter< RealType >( "SWD-pressure" );
- RealType eSWD = ( preSWD / ( rhoSWD * (gamma - 1) ) );
- //RealType eSWD = ( preSWD / (gamma - 1) ) + 0.5 * rhoSWD * pow(velSWDX,2)+pow(velSWDY,2);
+ //RealType eSWD = ( preSWD / ( rhoSWD * (gamma - 1) ) );
+ RealType eSWD = ( preSWD / (gamma - 1) ) + 0.5 * rhoSWD * (std::pow(velSWDX,2)+std::pow(velSWDY,2)+std::pow(velSWDZ,2));
RealType rhoNEU = parameters.getParameter< RealType >( "NEU-density" );
RealType velNEUX = parameters.getParameter< RealType >( "NEU-velocityX" );
RealType velNEUY = parameters.getParameter< RealType >( "NEU-velocityY" );
+ RealType velNEUZ = parameters.getParameter< RealType >( "NEU-velocityZ" );
RealType preNEU = parameters.getParameter< RealType >( "NEU-pressure" );
- RealType eNEU = ( preNEU / ( rhoNEU * (gamma - 1) ) );
- //RealType eNEU = ( preNEU / (gamma - 1) ) + 0.5 * rhoNEU * pow(velNEUX,2)+pow(velNEUY,2);
+ //RealType eNEU = ( preNEU / ( rhoNEU * (gamma - 1) ) );
+ RealType eNEU = ( preNEU / (gamma - 1) ) + 0.5 * rhoNEU * (std::pow(velNEUX,2)+std::pow(velNEUY,2)+std::pow(velNEUZ,2));
RealType rhoSEU = parameters.getParameter< RealType >( "SEU-density" );
RealType velSEUX = parameters.getParameter< RealType >( "SEU-velocityX" );
RealType velSEUY = parameters.getParameter< RealType >( "SEU-velocityY" );
+ RealType velSEUZ = parameters.getParameter< RealType >( "SEU-velocityZ" );
RealType preSEU = parameters.getParameter< RealType >( "SEU-pressure" );
- RealType eSEU = ( preSEU / ( rhoSEU * (gamma - 1) ) );
- //RealType eSEU = ( preSEU / (gamma - 1) ) + 0.5 * rhoSEU * pow(velSEUX,2)+pow(velSEUY,2);
+ //RealType eSEU = ( preSEU / ( rhoSEU * (gamma - 1) ) );
+ RealType eSEU = ( preSEU / (gamma - 1) ) + 0.5 * rhoSEU * (std::pow(velSEUX,2)+std::pow(velSEUY,2)+std::pow(velSEUZ,2));
RealType rhoNED = parameters.getParameter< RealType >( "NED-density" );
RealType velNEDX = parameters.getParameter< RealType >( "NED-velocityX" );
RealType velNEDY = parameters.getParameter< RealType >( "NED-velocityY" );
+ RealType velNEDZ = parameters.getParameter< RealType >( "NED-velocityZ" );
RealType preNED = parameters.getParameter< RealType >( "NED-pressure" );
- RealType eNED = ( preNED / ( rhoNED * (gamma - 1) ) );
- //RealType eNED = ( preNED / (gamma - 1) ) + 0.5 * rhoNED * pow(velNEDX,2)+pow(velNEDY,2);
+ //RealType eNED = ( preNED / ( rhoNED * (gamma - 1) ) );
+ RealType eNED = ( preNED / (gamma - 1) ) + 0.5 * rhoNED * (std::pow(velNEDX,2)+std::pow(velNEDY,2)+std::pow(velNEDZ,2));
RealType rhoSED = parameters.getParameter< RealType >( "SED-density" );
RealType velSEDX = parameters.getParameter< RealType >( "SED-velocityX" );
RealType velSEDY = parameters.getParameter< RealType >( "SED-velocityY" );
+ RealType velSEDZ = parameters.getParameter< RealType >( "SED-velocityZ" );
RealType preSED = parameters.getParameter< RealType >( "SED-pressure" );
- RealType eSED = ( preSED / ( rhoSED * (gamma - 1) ) );
- //RealType eSED = ( preSED / (gamma - 1) ) + 0.5 * rhoSED * pow(velSEDX,2)+pow(velSEDY,2);
+ //RealType eSED = ( preSED / ( rhoSED * (gamma - 1) ) );
+ RealType eSED = ( preSED / (gamma - 1) ) + 0.5 * rhoSED * (std::pow(velSEDX,2)+std::pow(velSEDY,2)+std::pow(velSEDZ,2));
RealType x0 = parameters.getParameter< RealType >( "riemann-border" );
int size = mesh->template getEntitiesCount< Cell >();
uRho->bind(mesh, dofs, 0);
@@ -166,7 +174,7 @@ setInitialCondition( const Config::ParameterContainer& parameters,
velocity->bind(mesh, data, size);
velocityX->bind(mesh, data, 2*size);
velocityY->bind(mesh, data, 3*size);
- velocityY->bind(mesh, data, 4*size);
+ velocityZ->bind(mesh, data, 4*size);
int count = std::pow(size, (1.0/3.0));
for(IndexType i = 0; i < count; i++)
for(IndexType j = 0; j < count; j++)
@@ -176,10 +184,12 @@ setInitialCondition( const Config::ParameterContainer& parameters,
(* uRho)[i*std::pow(count,2)+j*count+k] = rhoSWD;
(* uRhoVelocityX)[i*std::pow(count,2)+j*count+k] = rhoSWD * velSWDX;
(* uRhoVelocityY)[i*std::pow(count,2)+j*count+k] = rhoSWD * velSWDY;
+ (* uRhoVelocityZ)[i*std::pow(count,2)+j*count+k] = rhoSWD * velSWDZ;
(* uEnergy)[i*std::pow(count,2)+j*count+k] = eSWD;
- (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velSWDX,2)+std::pow(velSWDY,2));
+ (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velSWDX,2)+std::pow(velSWDY,2)+std::pow(velSWDZ,2));
(* velocityX)[i*std::pow(count,2)+j*count+k] = velSWDX;
(* velocityY)[i*std::pow(count,2)+j*count+k] = velSWDY;
+ (* velocityZ)[i*std::pow(count,2)+j*count+k] = velSWDZ;
(* pressure)[i*std::pow(count,2)+j*count+k] = preSWD;
}
else
@@ -188,10 +198,12 @@ setInitialCondition( const Config::ParameterContainer& parameters,
(* uRho)[i*std::pow(count,2)+j*count+k] = rhoSED;
(* uRhoVelocityX)[i*std::pow(count,2)+j*count+k] = rhoSED * velSEDX;
(* uRhoVelocityY)[i*std::pow(count,2)+j*count+k] = rhoSED * velSEDY;
+ (* uRhoVelocityZ)[i*std::pow(count,2)+j*count+k] = rhoSED * velSEDZ;
(* uEnergy)[i*std::pow(count,2)+j*count+k] = eSED;
- (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velSEDX,2)+std::pow(velSEDY,2));
+ (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velSEDX,2)+std::pow(velSEDY,2)+std::pow(velSEDZ,2));
(* velocityX)[i*std::pow(count,2)+j*count+k] = velSEDX;
(* velocityY)[i*std::pow(count,2)+j*count+k] = velSEDY;
+ (* velocityZ)[i*std::pow(count,2)+j*count+k] = velSEDZ;
(* pressure)[i*std::pow(count,2)+j*count+k] = preSED;
}
else
@@ -200,10 +212,12 @@ setInitialCondition( const Config::ParameterContainer& parameters,
(* uRho)[i*std::pow(count,2)+j*count+k] = rhoNWD;
(* uRhoVelocityX)[i*std::pow(count,2)+j*count+k] = rhoNWD * velNWDX;
(* uRhoVelocityY)[i*std::pow(count,2)+j*count+k] = rhoNWD * velNWDY;
+ (* uRhoVelocityZ)[i*std::pow(count,2)+j*count+k] = rhoNWD * velNWDZ;
(* uEnergy)[i*std::pow(count,2)+j*count+k] = eNWD;
- (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velNWDX,2)+std::pow(velNWDY,2));
+ (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velNWDX,2)+std::pow(velNWDY,2)+std::pow(velNWDZ,2));
(* velocityX)[i*std::pow(count,2)+j*count+k] = velNWDX;
(* velocityY)[i*std::pow(count,2)+j*count+k] = velNWDY;
+ (* velocityZ)[i*std::pow(count,2)+j*count+k] = velNWDZ;
(* pressure)[i*std::pow(count,2)+j*count+k] = preNWD;
}
else
@@ -212,10 +226,12 @@ setInitialCondition( const Config::ParameterContainer& parameters,
(* uRho)[i*std::pow(count,2)+j*count+k] = rhoNED;
(* uRhoVelocityX)[i*std::pow(count,2)+j*count+k] = rhoNED * velNEDX;
(* uRhoVelocityY)[i*std::pow(count,2)+j*count+k] = rhoNED * velNEDY;
+ (* uRhoVelocityZ)[i*std::pow(count,2)+j*count+k] = rhoNED * velNEDZ;
(* uEnergy)[i*std::pow(count,2)+j*count+k] = eNED;
- (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velNEDX,2)+std::pow(velNEDY,2));
+ (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velNEDX,2)+std::pow(velNEDY,2)+std::pow(velNEDZ,2));
(* velocityX)[i*std::pow(count,2)+j*count+k] = velNEDX;
(* velocityY)[i*std::pow(count,2)+j*count+k] = velNEDY;
+ (* velocityZ)[i*std::pow(count,2)+j*count+k] = velNEDZ;
(* pressure)[i*std::pow(count,2)+j*count+k] = preNED;
}
else
@@ -224,10 +240,12 @@ setInitialCondition( const Config::ParameterContainer& parameters,
(* uRho)[i*std::pow(count,2)+j*count+k] = rhoSWU;
(* uRhoVelocityX)[i*std::pow(count,2)+j*count+k] = rhoSWU * velSWUX;
(* uRhoVelocityY)[i*std::pow(count,2)+j*count+k] = rhoSWU * velSWUY;
+ (* uRhoVelocityZ)[i*std::pow(count,2)+j*count+k] = rhoSWU * velSWUZ;
(* uEnergy)[i*std::pow(count,2)+j*count+k] = eSWU;
- (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velSWUX,2)+std::pow(velSWUY,2));
+ (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velSWUX,2)+std::pow(velSWUY,2)+std::pow(velSWUZ,2));
(* velocityX)[i*std::pow(count,2)+j*count+k] = velSWUX;
(* velocityY)[i*std::pow(count,2)+j*count+k] = velSWUY;
+ (* velocityZ)[i*std::pow(count,2)+j*count+k] = velSWUZ;
(* pressure)[i*std::pow(count,2)+j*count+k] = preSWU;
}
else
@@ -236,10 +254,12 @@ setInitialCondition( const Config::ParameterContainer& parameters,
(* uRho)[i*std::pow(count,2)+j*count+k] = rhoSEU;
(* uRhoVelocityX)[i*std::pow(count,2)+j*count+k] = rhoSEU * velSEUX;
(* uRhoVelocityY)[i*std::pow(count,2)+j*count+k] = rhoSEU * velSEUY;
+ (* uRhoVelocityZ)[i*std::pow(count,2)+j*count+k] = rhoSEU * velSEUZ;
(* uEnergy)[i*std::pow(count,2)+j*count+k] = eSEU;
- (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velSEUX,2)+std::pow(velSEUY,2));
+ (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velSEUX,2)+std::pow(velSEUY,2)+std::pow(velSEUZ,2));
(* velocityX)[i*std::pow(count,2)+j*count+k] = velSEUX;
(* velocityY)[i*std::pow(count,2)+j*count+k] = velSEUY;
+ (* velocityZ)[i*std::pow(count,2)+j*count+k] = velSEUZ;
(* pressure)[i*std::pow(count,2)+j*count+k] = preSEU;
}
else
@@ -248,10 +268,12 @@ setInitialCondition( const Config::ParameterContainer& parameters,
(* uRho)[i*std::pow(count,2)+j*count+k] = rhoNWU;
(* uRhoVelocityX)[i*std::pow(count,2)+j*count+k] = rhoNWU * velNWUX;
(* uRhoVelocityY)[i*std::pow(count,2)+j*count+k] = rhoNWU * velNWUY;
+ (* uRhoVelocityZ)[i*std::pow(count,2)+j*count+k] = rhoNWU * velNWUZ;
(* uEnergy)[i*std::pow(count,2)+j*count+k] = eNWU;
- (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velNWUX,2)+std::pow(velNWUY,2));
+ (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velNWUX,2)+std::pow(velNWUY,2)+std::pow(velNWUZ,2));
(* velocityX)[i*std::pow(count,2)+j*count+k] = velNWUX;
(* velocityY)[i*std::pow(count,2)+j*count+k] = velNWUY;
+ (* velocityZ)[i*std::pow(count,2)+j*count+k] = velNWUZ;
(* pressure)[i*std::pow(count,2)+j*count+k] = preNWU;
}
else
@@ -260,10 +282,12 @@ setInitialCondition( const Config::ParameterContainer& parameters,
(* uRho)[i*std::pow(count,2)+j*count+k] = rhoNEU;
(* uRhoVelocityX)[i*std::pow(count,2)+j*count+k] = rhoNEU * velNEUX;
(* uRhoVelocityY)[i*std::pow(count,2)+j*count+k] = rhoNEU * velNEUY;
+ (* uRhoVelocityZ)[i*std::pow(count,2)+j*count+k] = rhoNEU * velNEUZ;
(* uEnergy)[i*std::pow(count,2)+j*count+k] = eNEU;
- (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velNEUX,2)+std::pow(velNEUY,2));
+ (* velocity)[i*std::pow(count,2)+j*count+k] = std::sqrt(std::pow(velNEUX,2)+std::pow(velNEUY,2)+std::pow(velNEUZ,2));
(* velocityX)[i*std::pow(count,2)+j*count+k] = velNEUX;
(* velocityY)[i*std::pow(count,2)+j*count+k] = velNEUY;
+ (* velocityZ)[i*std::pow(count,2)+j*count+k] = velNEUZ;
(* pressure)[i*std::pow(count,2)+j*count+k] = preNEU;
};
return true;
--
GitLab