diff --git a/src/GksGpu/BoundaryConditions/CreepingMassFlux.cu b/src/GksGpu/BoundaryConditions/CreepingMassFlux.cu
index 1e8179a8bc2e6b8f629eed3368f3e6fadee0dd3b..a164c73f5d1d8e1c1a9c7943acdd3ae485150778 100644
--- a/src/GksGpu/BoundaryConditions/CreepingMassFlux.cu
+++ b/src/GksGpu/BoundaryConditions/CreepingMassFlux.cu
@@ -136,7 +136,7 @@ bool CreepingMassFlux::isWall()
bool CreepingMassFlux::isFluxBC()
{
- return false;
+ return true;
}
bool CreepingMassFlux::secondCellsNeeded()
diff --git a/src/GksGpu/CellUpdate/CellUpdate.cu b/src/GksGpu/CellUpdate/CellUpdate.cu
index 71678138d3339ac8af9fe4077a9c711dc84ff90a..3b01f1d7b5bf310a8726cf8092000cd03093aa2c 100644
--- a/src/GksGpu/CellUpdate/CellUpdate.cu
+++ b/src/GksGpu/CellUpdate/CellUpdate.cu
@@ -82,6 +82,20 @@ __host__ __device__ inline void cellUpdateFunction(DataBaseStruct dataBase, Para
//PrimitiveVariables testPrim = toPrimitiveVariables(testCons, parameters.K);
//////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ //if( cellIndex == 415179 )
+ //{
+ // //printf( "rho = %14.4e | dRho = %14.4e \n", cons.rho , (one / cellVolume) * update.rho );
+ // //printf( "rhoU = %14.4e | dRhoU = %14.4e \n", cons.rhoU , (one / cellVolume) * update.rhoU );
+ // //printf( "rhoV = %14.4e | dRhoV = %14.4e \n", cons.rhoV , (one / cellVolume) * update.rhoV );
+ // //printf( "rhoW = %14.4e | dRhoW = %14.4e \n", cons.rhoW , (one / cellVolume) * update.rhoW );
+ // printf( "rhoE = %14.4e | dRhoE = %14.4e \n", cons.rhoE , (one / cellVolume) * update.rhoE );
+ // //printf( "rhoS1 = %14.4e | dRhoS1 = %14.4e \n", cons.rhoS_1, (one / cellVolume) * update.rhoS_1 );
+ // //printf( "rhoS2 = %14.4e | dRhoS2 = %14.4e \n", cons.rhoS_2, (one / cellVolume) * update.rhoS_2 );
+ // printf( "=================================================================\n" );
+ //}
+ //////////////////////////////////////////////////////////////////////////
+
cons = cons + (one / cellVolume) * update;
//////////////////////////////////////////////////////////////////////////
diff --git a/src/GksGpu/FluxComputation/AssembleFlux.cuh b/src/GksGpu/FluxComputation/AssembleFlux.cuh
index df7c21d58a8d6ece1c34e99325861e7d3f7011ce..28000cf81a89867bba081d19a95dc569863b426a 100644
--- a/src/GksGpu/FluxComputation/AssembleFlux.cuh
+++ b/src/GksGpu/FluxComputation/AssembleFlux.cuh
@@ -207,6 +207,8 @@ __host__ __device__ inline void computeTimeCoefficients(const PrimitiveVariables
{
real r = parameters.lambdaRef / facePrim.lambda;
+ //if( r < zero ) printf( "ERROR: %f/%f\n", parameters.lambdaRef, facePrim.lambda );
+
real mu;
if ( parameters.viscosityModel == ViscosityModel::constant ){
mu = parameters.mu;
@@ -539,19 +541,28 @@ __host__ __device__ inline void assembleFlux( const PrimitiveVariables& facePrim
flux_3.rhoS_2 = at[6] * momentU[0 + 1]
/ ( two * facePrim.lambda );
-
- real tauS = parameters.D * two * facePrim.lambda;
+ //////////////////////////////////////////////////////////////////////////
real dt = parameters.dt;
-
real timeCoefficientsScalar [3];
- timeCoefficientsScalar[0] = dt;
- timeCoefficientsScalar[1] = -tauS * dt;
- timeCoefficientsScalar[2] = c1o2 * dt * dt - tauS * dt;
+ {
+ real tauS = parameters.D1 * two * facePrim.lambda;
+ timeCoefficientsScalar[0] = dt;
+ timeCoefficientsScalar[1] = -tauS * dt;
+ timeCoefficientsScalar[2] = c1o2 * dt * dt - tauS * dt;
- flux.rhoS_1 = flux.rho * facePrim.S_1 + ( timeCoefficientsScalar[1] * flux_2.rhoS_1 + timeCoefficientsScalar[2] * flux_3.rhoS_1 ) * parameters.dx * parameters.dx * facePrim.rho;
- flux.rhoS_2 = flux.rho * facePrim.S_2 + ( timeCoefficientsScalar[1] * flux_2.rhoS_2 + timeCoefficientsScalar[2] * flux_3.rhoS_2 ) * parameters.dx * parameters.dx * facePrim.rho;
+ flux.rhoS_1 = flux.rho * facePrim.S_1 + ( timeCoefficientsScalar[1] * flux_2.rhoS_1 + timeCoefficientsScalar[2] * flux_3.rhoS_1 ) * parameters.dx * parameters.dx * facePrim.rho;
+ }
+
+ {
+ real tauS = parameters.D2 * two * facePrim.lambda;
+ timeCoefficientsScalar[0] = dt;
+ timeCoefficientsScalar[1] = -tauS * dt;
+ timeCoefficientsScalar[2] = c1o2 * dt * dt - tauS * dt;
+
+ flux.rhoS_2 = flux.rho * facePrim.S_2 + ( timeCoefficientsScalar[1] * flux_2.rhoS_2 + timeCoefficientsScalar[2] * flux_3.rhoS_2 ) * parameters.dx * parameters.dx * facePrim.rho;
+ }
#endif // USE_PASSIVE_SCALAR
}
diff --git a/src/GksGpu/FluxComputation/FluxComputation.cu b/src/GksGpu/FluxComputation/FluxComputation.cu
index a397b5af5fb47a79b111f04eb262483de5fe5619..f9344fe2737e9dd0a481b72e554aa628c5e9c234 100644
--- a/src/GksGpu/FluxComputation/FluxComputation.cu
+++ b/src/GksGpu/FluxComputation/FluxComputation.cu
@@ -114,6 +114,9 @@ __host__ __device__ inline void fluxFunction(DataBaseStruct dataBase, Parameters
direction = dataBase.faceOrientation[ faceIndex ];
+ parameters.D1 = parameters.D;
+ parameters.D2 = parameters.D;
+
//////////////////////////////////////////////////////////////////////////
if( parameters.useSpongeLayer )
@@ -259,17 +262,19 @@ __host__ __device__ inline void fluxFunction(DataBaseStruct dataBase, Parameters
//////////////////////////////////////////////////////////////////////////
+ parameters.D1 = parameters.D;
+ parameters.D2 = parameters.D;
+
if(parameters.usePassiveScalarLimiter){
#ifdef USE_PASSIVE_SCALAR
- if( facePrim.S_1 < zero )
- {
- parameters.D += - parameters.passiveScalarLimiter * facePrim.S_1;
- }
- if( facePrim.S_1 > one )
- {
- parameters.D += parameters.passiveScalarLimiter * ( facePrim.S_1 - one );
- }
+ if( facePrim.S_1 < zero ) parameters.D1 += - parameters.passiveScalarLimiter * facePrim.S_1;
+ if( facePrim.S_1 > one ) parameters.D1 += parameters.passiveScalarLimiter * ( facePrim.S_1 - one );
+
+ parameters.D2 = parameters.D1;
+
+ if( facePrim.S_2 < zero ) parameters.D2 += - real(0.1)*parameters.passiveScalarLimiter * facePrim.S_2;
+ if( facePrim.S_2 > one ) parameters.D2 += real(0.1)*parameters.passiveScalarLimiter * ( facePrim.S_2 - one );
#endif // USE_PASSIVE_SCALAR
}
diff --git a/src/GksGpu/Parameters/Parameters.h b/src/GksGpu/Parameters/Parameters.h
index db640b877d2bd605392a218852d01689146c3118..2ac77419aabaaafef836f647cd8ceffd9cdfe04c 100644
--- a/src/GksGpu/Parameters/Parameters.h
+++ b/src/GksGpu/Parameters/Parameters.h
@@ -18,6 +18,8 @@ struct VF_PUBLIC Parameters
real K = real(2.0);
real Pr = real(1.0);
real D = real(0.01);
+ real D1 = real(0.01);
+ real D2 = real(0.01);
real dt = real(0.01);
real dx = real(0.01);
diff --git a/targets/apps/GKS/Flame7cm/Flame7cm.cpp b/targets/apps/GKS/Flame7cm/Flame7cm.cpp
index c87917d2a8c9702426d62a86c1ab1df291a27f5b..490ab1b14d01cbcb59ff7649e4b01f955a46aedf 100644
--- a/targets/apps/GKS/Flame7cm/Flame7cm.cpp
+++ b/targets/apps/GKS/Flame7cm/Flame7cm.cpp
@@ -63,7 +63,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
{
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- uint nx = 128;
+ uint nx = 256;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -90,7 +90,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
real cs = sqrt( ( ( K + 5.0 ) / ( K + 3.0 ) ) / ( 2.0 * prim.lambda ) );
- real CFL = 0.125;
+ real CFL = 0.06125;
real dt = CFL * ( dx / ( ( U + cs ) * ( one + ( two * mu ) / ( U * dx * rho ) ) ) );
@@ -138,7 +138,7 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
parameters.useReactionLimiter = true;
parameters.useTemperatureLimiter = true;
parameters.usePassiveScalarLimiter = true;
- parameters.useSmagorinsky = false;
+ parameters.useSmagorinsky = true;
parameters.reactionLimiter = 1.005;
@@ -380,17 +380,17 @@ void thermalCavity( std::string path, std::string simulationName, uint restartIt
}
if(
- //( iter >= 144820 && iter % 1 == 0 ) ||
+ //( iter >= 39360 && iter % 1 == 0 ) ||
( iter % 10000 == 0 )
)
{
dataBase->copyDataDeviceToHost();
-
writeVtkXML( dataBase, parameters, 0, path + simulationName + "_" + std::to_string( iter ) );
}
- if( iter % 10000 == 0 )
+ if( iter % 10000 == 0 /*|| iter == 39000*/)
{
+ dataBase->copyDataDeviceToHost();
Restart::writeRestart( dataBase, path + simulationName + "_" + std::to_string( iter ), iter );
}
@@ -442,7 +442,7 @@ int main( int argc, char* argv[])
try
{
uint restartIter = INVALID_INDEX;
- //uint restartIter = 140000;
+ //uint restartIter = 30000;
if( argc > 1 ) restartIter = atoi( argv[1] );