diff --git a/src/GksGpu/FluxComputation/FluxComputation.cu b/src/GksGpu/FluxComputation/FluxComputation.cu
index c3692465bbd8b7a7bfc7a806982f79cf2292bd5c..dd6ae0f5e9e84e19091011e48cc460110bd5ede6 100644
--- a/src/GksGpu/FluxComputation/FluxComputation.cu
+++ b/src/GksGpu/FluxComputation/FluxComputation.cu
@@ -183,10 +183,84 @@ __host__ __device__ inline void fluxFunction(DataBaseStruct dataBase, Parameters
//parameters.mu += getTurbulentViscositySmagorinsky( parameters, facePrim, gradN, gradT1, gradT2 );
//parameters.D = parameters.mu;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ {
+ real k = parameters.mu / parameters.Pr;
+
+ real dUdx1 = ( gradN.rhoU - facePrim.U * gradN.rho );
+ real dUdx2 = ( gradT1.rhoU - facePrim.U * gradT1.rho );
+ real dUdx3 = ( gradT2.rhoU - facePrim.U * gradT2.rho );
+
+ real dVdx1 = ( gradN.rhoV - facePrim.V * gradN.rho );
+ real dVdx2 = ( gradT1.rhoV - facePrim.V * gradT1.rho );
+ real dVdx3 = ( gradT2.rhoV - facePrim.V * gradT2.rho );
+
+ real dWdx1 = ( gradN.rhoW - facePrim.W * gradN.rho );
+ real dWdx2 = ( gradT1.rhoW - facePrim.W * gradT1.rho );
+ real dWdx3 = ( gradT2.rhoW - facePrim.W * gradT2.rho );
+
+ real dEdx1 = ( gradN.rhoE - facePrim.W * gradN.rho );
+ real dEdx2 = ( gradT1.rhoE - facePrim.W * gradT1.rho );
+ real dEdx3 = ( gradT2.rhoE - facePrim.W * gradT2.rho );
+
+ real dTdx1 = dEdx1 - two * facePrim.U * dUdx1 - two * facePrim.V * dVdx1 - two * facePrim.W * dWdx1;
+ real dTdx2 = dEdx2 - two * facePrim.U * dUdx2 - two * facePrim.V * dVdx2 - two * facePrim.W * dWdx2;
+ real dTdx3 = dEdx3 - two * facePrim.U * dUdx3 - two * facePrim.V * dVdx3 - two * facePrim.W * dWdx3;
+
+ // this one works for some time
+ real S = parameters.dx * parameters.dx * ( fabsf(dTdx1) + fabsf(dTdx2) + fabsf(dTdx3) );
+ k += real(0.00002) / real(0.015625) * S;
+
+
+ //real S = parameters.dx * ( fabsf(dTdx1) + fabsf(dTdx2) + fabsf(dTdx3) );
+ //k += real(0.00001) * real(0.0025) * S * S;
+
+ //real S = parameters.dx * parameters.dx * ( dTdx1 * dTdx1 + dTdx2 * dTdx2 + dTdx3 * dTdx3 );
+ //k += real(0.00001) * real(0.001) * S;
+
+ parameters.Pr = parameters.mu / k;
+
+ //if( parameters.Pr < 0.1 ) printf("%f\n", S);
+ }
}
//////////////////////////////////////////////////////////////////////////
+
+ //{
+ //#ifdef USE_PASSIVE_SCALAR
+ // if( facePrim.S_1 < zero )
+ // {
+ // parameters.D += - real(0.1) * facePrim.S_1;
+ // }
+ // if( facePrim.S_1 > one )
+ // {
+ // parameters.D += real(0.1) * ( facePrim.S_1 - one );
+ // }
+
+ //#endif // USE_PASSIVE_SCALAR
+ //}
+
+ //////////////////////////////////////////////////////////////////////////
+
+ //{
+ //#ifdef USE_PASSIVE_SCALAR
+ // if( facePrim.S_1 < zero )
+ // {
+ // parameters.D += - real(0.1) * facePrim.S_1;
+ // }
+ // if( facePrim.S_1 > one )
+ // {
+ // parameters.D += real(0.1) * ( facePrim.S_1 - one );
+ // }
+
+ //#endif // USE_PASSIVE_SCALAR
+ //}
+
+ //////////////////////////////////////////////////////////////////////////
+
{
ConservedVariables flux;
diff --git a/src/GksGpu/FluxComputation/Smagorinsky.cuh b/src/GksGpu/FluxComputation/Smagorinsky.cuh
index 99ebab67684b9ce4e4dcab06f54602d039b7635a..53b924c347c2d532cc2b557f71d21399202e10fa 100644
--- a/src/GksGpu/FluxComputation/Smagorinsky.cuh
+++ b/src/GksGpu/FluxComputation/Smagorinsky.cuh
@@ -21,17 +21,15 @@ inline __host__ __device__ real getTurbulentViscositySmagorinsky(const Parameter
{
// See FDS 6 Technical Reference Guide, Section 4.2.8
- real dUdx1 = ( gradX1.rhoU - facePrim.U * gradX1.rho ) / facePrim.rho;
- real dUdx2 = ( gradX2.rhoU - facePrim.U * gradX2.rho ) / facePrim.rho;
- real dUdx3 = ( gradX3.rhoU - facePrim.U * gradX3.rho ) / facePrim.rho;
-
- real dVdx1 = ( gradX1.rhoV - facePrim.V * gradX1.rho ) / facePrim.rho;
- real dVdx2 = ( gradX2.rhoV - facePrim.V * gradX2.rho ) / facePrim.rho;
- real dVdx3 = ( gradX3.rhoV - facePrim.V * gradX3.rho ) / facePrim.rho;
-
- real dWdx1 = ( gradX1.rhoW - facePrim.W * gradX1.rho ) / facePrim.rho;
- real dWdx2 = ( gradX2.rhoW - facePrim.W * gradX2.rho ) / facePrim.rho;
- real dWdx3 = ( gradX3.rhoW - facePrim.W * gradX3.rho ) / facePrim.rho;
+ real dUdx1 = ( gradX1.rhoU - facePrim.U * gradX1.rho )/* / facePrim.rho*/;
+ real dUdx2 = ( gradX2.rhoU - facePrim.U * gradX2.rho )/* / facePrim.rho*/;
+ real dUdx3 = ( gradX3.rhoU - facePrim.U * gradX3.rho )/* / facePrim.rho*/;
+ real dVdx1 = ( gradX1.rhoV - facePrim.V * gradX1.rho )/* / facePrim.rho*/;
+ real dVdx2 = ( gradX2.rhoV - facePrim.V * gradX2.rho )/* / facePrim.rho*/;
+ real dVdx3 = ( gradX3.rhoV - facePrim.V * gradX3.rho )/* / facePrim.rho*/;
+ real dWdx1 = ( gradX1.rhoW - facePrim.W * gradX1.rho )/* / facePrim.rho*/;
+ real dWdx2 = ( gradX2.rhoW - facePrim.W * gradX2.rho )/* / facePrim.rho*/;
+ real dWdx3 = ( gradX3.rhoW - facePrim.W * gradX3.rho )/* / facePrim.rho*/;
real S11sq = dUdx1*dUdx1;
real S22sq = dVdx2*dVdx2;
diff --git a/targets/apps/GKS/PoolFire/PoolFire.cpp b/targets/apps/GKS/PoolFire/PoolFire.cpp
index df454e0093d7afa116adcb6d00bb01df196e6d27..00c77a1e3a6420695c358a05e6195294ab69e825 100644
--- a/targets/apps/GKS/PoolFire/PoolFire.cpp
+++ b/targets/apps/GKS/PoolFire/PoolFire.cpp
@@ -99,6 +99,8 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
*logging::out << logging::Logger::INFO_HIGH << "cs = " << cs << " m/s\n";
*logging::out << logging::Logger::INFO_HIGH << "mu = " << mu << " kg/sm\n";
+ *logging::out << logging::Logger::INFO_HIGH << "F_rho = " << U * rho * dt * 1000 << " kg/m^3\n";
+
//////////////////////////////////////////////////////////////////////////
Parameters parameters;
@@ -215,10 +217,10 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
real openBoundaryVelocityLimiter = 1.0;
- //SPtr<BoundaryCondition> bcMX = std::make_shared<Open>( dataBase, prim, openBoundaryVelocityLimiter );
- //SPtr<BoundaryCondition> bcPX = std::make_shared<Open>( dataBase, prim, openBoundaryVelocityLimiter );
- SPtr<BoundaryCondition> bcMX = std::make_shared<AdiabaticWall>( dataBase, Vec3(0, 0, 0), true );
- SPtr<BoundaryCondition> bcPX = std::make_shared<AdiabaticWall>( dataBase, Vec3(0, 0, 0), true );
+ SPtr<BoundaryCondition> bcMX = std::make_shared<Open>( dataBase, prim, openBoundaryVelocityLimiter );
+ SPtr<BoundaryCondition> bcPX = std::make_shared<Open>( dataBase, prim, openBoundaryVelocityLimiter );
+ //SPtr<BoundaryCondition> bcMX = std::make_shared<AdiabaticWall>( dataBase, Vec3(0, 0, 0), true );
+ //SPtr<BoundaryCondition> bcPX = std::make_shared<AdiabaticWall>( dataBase, Vec3(0, 0, 0), true );
//SPtr<BoundaryCondition> bcMX = std::make_shared<MassCompensation>( dataBase, rho, U, prim.lambda );
//SPtr<BoundaryCondition> bcPX = std::make_shared<MassCompensation>( dataBase, rho, U, prim.lambda );
//SPtr<BoundaryCondition> bcMX = std::make_shared<IsothermalWall>( dataBase, Vec3(0, 0, 0), prim.lambda, false );
@@ -272,10 +274,10 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
//SPtr<BoundaryCondition> bcPZ = std::make_shared<Open>( dataBase, prim );
//SPtr<BoundaryCondition> bcPZ = std::make_shared<Extrapolation>( dataBase );
- SPtr<BoundaryCondition> bcPZ = std::make_shared<AdiabaticWall>( dataBase, Vec3(0, 0, 0), true );
- //SPtr<BoundaryCondition> bcPZ = std::make_shared<Pressure2>( dataBase, c1o2 * prim.rho / prim.lambda );
+ //SPtr<BoundaryCondition> bcPZ = std::make_shared<AdiabaticWall>( dataBase, Vec3(0, 0, 0), true );
+ SPtr<BoundaryCondition> bcPZ = std::make_shared<Pressure2>( dataBase, c1o2 * prim.rho / prim.lambda );
- bcMZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z < 0.0 || ( std::sqrt(center.x*center.x + center.y*center.y) < 0.6 && center.z < 0.1); } );
+ bcMZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z < 0.0; } );
bcPZ->findBoundaryCells( meshAdapter, true, [&](Vec3 center){ return center.z > H ; } );
//////////////////////////////////////////////////////////////////////////
@@ -287,7 +289,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
burner->findBoundaryCells( meshAdapter, false, [&](Vec3 center){
- return center.z < 0.0 && std::sqrt(center.x*center.x + center.y*center.y) < 0.5;
+ return center.z < 0.0 && std::sqrt(center.x*center.x) < 0.5 && std::sqrt(center.y*center.y) < 0.5 * dx;
} );
//////////////////////////////////////////////////////////////////////////
@@ -368,7 +370,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
cupsAnalyzer.start();
- for( uint iter = startIter + 1; iter <= 20010; iter++ )
+ for( uint iter = startIter + 1; iter <= 2000000; iter++ )
{
uint runUpTime = 10000;
@@ -386,11 +388,11 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
//parameters.dt = 0.2 * dt + ( dt - 0.2 * dt ) * ( real(iter) / 40000.0 );
}
- if( iter == 5001 )
- {
- parameters.enableReaction = false;
- std::dynamic_pointer_cast<CreepingMassFlux>(burner)->velocity = -1.0;
- }
+ //if( iter == 5001 )
+ //{
+ // parameters.enableReaction = false;
+ // std::dynamic_pointer_cast<CreepingMassFlux>(burner)->velocity = -1.0;
+ //}
cupsAnalyzer.run( iter );
@@ -399,7 +401,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
TimeStepping::nestedTimeStep(dataBase, parameters, 0);
if(
- ( iter >= 20000 && iter % 1 == 0 ) ||
+ //( iter >= 20000 && iter % 1 == 0 ) ||
( iter % 1000 == 0 )
)
{
@@ -408,7 +410,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
writeVtkXML( dataBase, parameters, 0, path + simulationName + "_" + std::to_string( iter ) );
}
- if( iter % 10000 == 0 )
+ if( iter % 1000 == 0 )
{
Restart::writeRestart( dataBase, path + simulationName + "_" + std::to_string( iter ), iter );
}
@@ -450,7 +452,7 @@ int main( int argc, char* argv[])
try
{
uint restartIter = INVALID_INDEX;
- //uint restartIter = 20000;
+ //uint restartIter = 35000;
if( argc > 1 ) restartIter = atoi( argv[1] );