Skip to content
Snippets Groups Projects
Commit e905aadb authored by Henrik Asmuth's avatar Henrik Asmuth
Browse files

Dummy Stress BC with samplingIndices now working

parent e27eaa88
No related branches found
No related tags found
1 merge request!98Merge new StressBC, new Probes and various small changes from UU
Hide whitespace changes
Inline Side-by-side
Showing
with 259 additions and 441 deletions
apps/gpu/LBM/BoundaryLayer/BoundaryLayer.cpp
+ 3
− 3
View file @ e905aadb
......@@ -185,10 +185,10 @@ void multipleLevel(const std::string& configPath)
uint samplingOffset = 1;
// gridBuilder->setVelocityBoundaryCondition(SideType::PZ, 0.0, 0.0, 0.0);
// gridBuilder->setVelocityBoundaryCondition(SideType::MZ, 0.0, 0.0, 0.0);
gridBuilder->setStressBoundaryCondition(SideType::PZ, 0.0, 0.0, 0.0, samplingOffset);
// gridBuilder->setVelocityBoundaryCondition(SideType::PZ, 0.05, 0.0, 0.0);
gridBuilder->setStressBoundaryCondition(SideType::MZ, 0.0, 0.0, 1.0, samplingOffset);
gridBuilder->setVelocityBoundaryCondition(SideType::PZ, 0.05, 0.0, 0.0);
gridBuilder->setSlipBoundaryCondition(SideType::MZ, 0.0, 0.0, 0.0);
// gridBuilder->setSlipBoundaryCondition(SideType::MZ, 0.0, 0.0, 0.0);
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
SPtr<CudaMemoryManager> cudaMemoryManager = CudaMemoryManager::make(para);
......
src/gpu/GridGenerator/grid/BoundaryConditions/BoundaryCondition.cpp
+ 23
− 0
View file @ e905aadb
......@@ -42,6 +42,8 @@ bool gg::BoundaryCondition::isSide( SideType side ) const
return this->side->whoAmI() == side;
}
//////////////////////////////////////////////////////////////////////////
void VelocityBoundaryCondition::setVelocityProfile(
SPtr<Grid> grid, std::function<void(real, real, real, real &, real &, real &)> velocityProfile)
{
......@@ -55,6 +57,8 @@ void VelocityBoundaryCondition::setVelocityProfile(
}
}
//////////////////////////////////////////////////////////////////////////
void GeometryBoundaryCondition::setTangentialVelocityForPatch(SPtr<Grid> grid, uint patch,
real p1x, real p1y, real p1z,
real p2x, real p2y, real p2z,
......@@ -103,3 +107,22 @@ void GeometryBoundaryCondition::setTangentialVelocityForPatch(SPtr<Grid> grid, u
}
}
//////////////////////////////////////////////////////////////////////////
void StressBoundaryCondition::fillSamplingIndices(std::vector<SPtr<Grid> > grid, uint level, uint samplingOffset)
{
for( uint i = 0; i < this->indices.size(); i++ )
{
real x, y, z;
grid[level]->transIndexToCoords(this->indices[i], x, y, z);
real x_sampling = x + this->getNormalx(i)*samplingOffset*grid[level]->getDelta();
real y_sampling = y + this->getNormaly(i)*samplingOffset*grid[level]->getDelta();
real z_sampling = z + this->getNormalz(i)*samplingOffset*grid[level]->getDelta();
this->velocitySamplingIndices.push_back( grid[level]->transCoordToIndex(x_sampling, y_sampling, z_sampling) );
}
}
src/gpu/GridGenerator/grid/BoundaryConditions/BoundaryCondition.h
+ 2
− 0
View file @ e905aadb
......@@ -177,6 +177,8 @@ public:
real getNormalx(uint index) { return this->normalXList[index]; }
real getNormaly(uint index) { return this->normalYList[index]; }
real getNormalz(uint index) { return this->normalZList[index]; }
void fillSamplingIndices(std::vector<SPtr<Grid> > grid, uint level, uint samplingOffset);
};
//////////////////////////////////////////////////////////////////////////
......
src/gpu/GridGenerator/grid/BoundaryConditions/Side.cpp
+ 3
− 3
View file @ e905aadb
......@@ -106,11 +106,11 @@ void Side::setStressSamplingIndices(SPtr<BoundaryCondition> boundaryCondition, S
real nz = z;
if (boundaryCondition->side->getCoordinate() == X_INDEX)
nx = boundaryCondition->side->getDirection() * stressBoundaryCondition->samplingOffset * grid->getDelta() + x;
nx = -boundaryCondition->side->getDirection() * stressBoundaryCondition->samplingOffset * grid->getDelta() + x;
if (boundaryCondition->side->getCoordinate() == Y_INDEX)
ny = boundaryCondition->side->getDirection() * stressBoundaryCondition->samplingOffset * grid->getDelta() + y;
ny = -boundaryCondition->side->getDirection() * stressBoundaryCondition->samplingOffset * grid->getDelta() + y;
if (boundaryCondition->side->getCoordinate() == Z_INDEX)
nz = boundaryCondition->side->getDirection() * stressBoundaryCondition->samplingOffset * grid->getDelta() + z;
nz = -boundaryCondition->side->getDirection() * stressBoundaryCondition->samplingOffset * grid->getDelta() + z;
uint samplingIndex = grid->transCoordToIndex(nx, ny, nz);
stressBoundaryCondition->velocitySamplingIndices.push_back(samplingIndex);
......
src/gpu/GridGenerator/grid/GridBuilder/LevelGridBuilder.cpp
+ 3
− 1
View file @ e905aadb
......@@ -99,6 +99,7 @@ void LevelGridBuilder::setStressBoundaryCondition(SideType sideType, real nomalX
stressBoundaryCondition->side->addIndices(grids, 0, stressBoundaryCondition);
stressBoundaryCondition->fillStressNormalLists();
// stressBoundaryCondition->fillSamplingIndices(grids, 0, samplingOffset); //redundant with Side::setStressSamplingIndices but potentially a better approach for cases with complex geometries
boundaryConditions[0]->stressBoundaryConditions.push_back(stressBoundaryCondition);
......@@ -394,7 +395,8 @@ void LevelGridBuilder::getStressValues(real* normalX, real* normalY, real* norma
{
indices[allIndicesCounter] = grids[level]->getSparseIndex(boundaryCondition->indices[index]) + 1;
samplingIndices[allIndicesCounter] = grids[level]->getSparseIndex(boundaryCondition->velocitySamplingIndices[index]) + 1;
if(index==100)
printf("IDX %u \t IDXs %u\n",boundaryCondition->indices[index], boundaryCondition->velocitySamplingIndices[index] );
normalX[allIndicesCounter] = boundaryCondition->getNormalx(index);
normalY[allIndicesCounter] = boundaryCondition->getNormaly(index);
normalZ[allIndicesCounter] = boundaryCondition->getNormalz(index);
......
src/gpu/VirtualFluids_GPU/GPU/StressBCs27.cu
+ 225
− 434
View file @ e905aadb
......@@ -276,7 +276,7 @@ extern "C" __global__ void QStressDeviceComp27(real* DD,
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//Compute incoming f's with zero wall velocity
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
real VeloX=0.05, VeloY=0.0, VeloZ=0.0;
// real VeloX=0.0, VeloY=0.0, VeloZ=0.0;
// if(k==3071){ // 3071 9120 9215
// printf("======================================\n");
......@@ -340,621 +340,412 @@ extern "C" __global__ void QStressDeviceComp27(real* DD,
// printf("TSW \t %f \n\n", f_TSW);
// }
//ToDo anders Klammern
// incoming f's from bounce back
real f_E_in, f_W_in, f_N_in, f_S_in, f_T_in, f_B_in, f_NE_in, f_SW_in, f_SE_in, f_NW_in, f_TE_in, f_BW_in, f_BE_in,
f_TW_in, f_TN_in, f_BS_in, f_BN_in, f_TS_in, f_TNE_in, f_TSW_in, f_TSE_in, f_TNW_in, f_BNE_in, f_BSW_in, f_BSE_in, f_BNW_in;
// momentum exchanged with wall at rest
real wallMomentumX = 0.0, wallMomentumY = 0.0, wallMomentumZ = 0.0;
q = q_dirE[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c2o27* (drho/*+three*( vx1 )*/+c9o2*( vx1 )*( vx1 ) * (c1o1 + drho)-cu_sq);
(D.f[dirW])[kw]=(c1o1-q)/(c1o1+q)*(f_E-f_W+(f_E+f_W-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_E+f_W)-c6o1*c2o27*( VeloX ))/(c1o1+q) - c2o27 * drho;
//(D.f[dirW])[kw]=zero;
f_W_in=(c1o1-q)/(c1o1+q)*(f_E-f_W+(f_E+f_W-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_E+f_W))/(c1o1+q) - c2o27 * drho;
wallMomentumX += f_E+f_W_in;
}
q = q_dirW[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c2o27* (drho/*+three*(-vx1 )*/+c9o2*(-vx1 )*(-vx1 ) * (c1o1 + drho)-cu_sq);
(D.f[dirE])[ke]=(c1o1-q)/(c1o1+q)*(f_W-f_E+(f_W+f_E-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_W+f_E)-c6o1*c2o27*(-VeloX ))/(c1o1+q) - c2o27 * drho;
//(D.f[dirE])[ke]=zero;
f_E_in=(c1o1-q)/(c1o1+q)*(f_W-f_E+(f_W+f_E-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_W+f_E))/(c1o1+q) - c2o27 * drho;
wallMomentumX -= f_W+f_E_in;
}
q = q_dirN[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c2o27* (drho/*+three*( vx2 )*/+c9o2*( vx2 )*( vx2 ) * (c1o1 + drho)-cu_sq);
(D.f[dirS])[ks]=(c1o1-q)/(c1o1+q)*(f_N-f_S+(f_N+f_S-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_N+f_S)-c6o1*c2o27*( VeloY ))/(c1o1+q) - c2o27 * drho;
//(D.f[dirS])[ks]=zero;
f_S_in=(c1o1-q)/(c1o1+q)*(f_N-f_S+(f_N+f_S-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_N+f_S))/(c1o1+q) - c2o27 * drho;
wallMomentumY += f_N+f_S_in;
}
q = q_dirS[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c2o27* (drho/*+three*( -vx2 )*/+c9o2*( -vx2 )*( -vx2 ) * (c1o1 + drho)-cu_sq);
(D.f[dirN])[kn]=(c1o1-q)/(c1o1+q)*(f_S-f_N+(f_S+f_N-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_S+f_N)-c6o1*c2o27*(-VeloY ))/(c1o1+q) - c2o27 * drho;
//(D.f[dirN])[kn]=zero;
f_N_in=(c1o1-q)/(c1o1+q)*(f_S-f_N+(f_S+f_N-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_S+f_N))/(c1o1+q) - c2o27 * drho;
wallMomentumY -= f_S+f_N_in;
}
q = q_dirT[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c2o27* (drho/*+three*( vx3)*/+c9o2*( vx3)*( vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirB])[kb]=(c1o1-q)/(c1o1+q)*(f_T-f_B+(f_T+f_B-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_T+f_B)-c6o1*c2o27*( VeloZ ))/(c1o1+q) - c2o27 * drho;
//(D.f[dirB])[kb]=one;
f_B_in=(c1o1-q)/(c1o1+q)*(f_T-f_B+(f_T+f_B-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_T+f_B))/(c1o1+q) - c2o27 * drho;
wallMomentumZ += f_T+f_B_in;
}
q = q_dirB[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c2o27* (drho/*+three*( -vx3)*/+c9o2*( -vx3)*( -vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirT])[kt]=(c1o1-q)/(c1o1+q)*(f_B-f_T+(f_B+f_T-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_B+f_T)-c6o1*c2o27*(-VeloZ ))/(c1o1+q) - c2o27 * drho;
//(D.f[dirT])[kt]=zero;
f_T_in=(c1o1-q)/(c1o1+q)*(f_B-f_T+(f_B+f_T-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_B+f_T))/(c1o1+q) - c2o27 * drho;
wallMomentumZ -= f_B+f_T_in;
}
q = q_dirNE[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o54* (drho/*+three*( vx1+vx2 )*/+c9o2*( vx1+vx2 )*( vx1+vx2 ) * (c1o1 + drho)-cu_sq);
(D.f[dirSW])[ksw]=(c1o1-q)/(c1o1+q)*(f_NE-f_SW+(f_NE+f_SW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_NE+f_SW)-c6o1*c1o54*(VeloX+VeloY))/(c1o1+q) - c1o54 * drho;
//(D.f[dirSW])[ksw]=zero;
f_SW_in=(c1o1-q)/(c1o1+q)*(f_NE-f_SW+(f_NE+f_SW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_NE+f_SW))/(c1o1+q) - c1o54 * drho;
wallMomentumX += f_NE+f_SW_in;
wallMomentumY += f_NE+f_SW_in;
}
q = q_dirSW[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o54* (drho/*+three*(-vx1-vx2 )*/+c9o2*(-vx1-vx2 )*(-vx1-vx2 ) * (c1o1 + drho)-cu_sq);
(D.f[dirNE])[kne]=(c1o1-q)/(c1o1+q)*(f_SW-f_NE+(f_SW+f_NE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_SW+f_NE)-c6o1*c1o54*(-VeloX-VeloY))/(c1o1+q) - c1o54 * drho;
//(D.f[dirNE])[kne]=zero;
f_NE_in=(c1o1-q)/(c1o1+q)*(f_SW-f_NE+(f_SW+f_NE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_SW+f_NE))/(c1o1+q) - c1o54 * drho;
wallMomentumX -= f_SW+f_NE_in;
wallMomentumY -= f_SW+f_NE_in;
}
q = q_dirSE[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o54* (drho/*+three*( vx1-vx2 )*/+c9o2*( vx1-vx2 )*( vx1-vx2 ) * (c1o1 + drho)-cu_sq);
(D.f[dirNW])[knw]=(c1o1-q)/(c1o1+q)*(f_SE-f_NW+(f_SE+f_NW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_SE+f_NW)-c6o1*c1o54*( VeloX-VeloY))/(c1o1+q) - c1o54 * drho;
//(D.f[dirNW])[knw]=zero;
f_NW_in=(c1o1-q)/(c1o1+q)*(f_SE-f_NW+(f_SE+f_NW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_SE+f_NW))/(c1o1+q) - c1o54 * drho;
wallMomentumX += f_SE+f_NW_in;
wallMomentumY -= f_SE+f_NW_in;
}
q = q_dirNW[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o54* (drho/*+three*(-vx1+vx2 )*/+c9o2*(-vx1+vx2 )*(-vx1+vx2 ) * (c1o1 + drho)-cu_sq);
(D.f[dirSE])[kse]=(c1o1-q)/(c1o1+q)*(f_NW-f_SE+(f_NW+f_SE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_NW+f_SE)-c6o1*c1o54*(-VeloX+VeloY))/(c1o1+q) - c1o54 * drho;
//(D.f[dirSE])[kse]=zero;
f_SE_in=(c1o1-q)/(c1o1+q)*(f_NW-f_SE+(f_NW+f_SE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_NW+f_SE))/(c1o1+q) - c1o54 * drho;
wallMomentumX -= f_NW+f_SE_in;
wallMomentumY += f_NW+f_SE_in;
}
q = q_dirTE[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o54* (drho/*+three*( vx1 +vx3)*/+c9o2*( vx1 +vx3)*( vx1 +vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirBW])[kbw]=(c1o1-q)/(c1o1+q)*(f_TE-f_BW+(f_TE+f_BW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TE+f_BW)-c6o1*c1o54*( VeloX+VeloZ))/(c1o1+q) - c1o54 * drho;
//(D.f[dirBW])[kbw]=zero;
f_BW_in=(c1o1-q)/(c1o1+q)*(f_TE-f_BW+(f_TE+f_BW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TE+f_BW))/(c1o1+q) - c1o54 * drho;
wallMomentumX += f_TE+f_BW_in;
wallMomentumZ += f_TE+f_BW_in;
}
q = q_dirBW[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o54* (drho/*+three*(-vx1 -vx3)*/+c9o2*(-vx1 -vx3)*(-vx1 -vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirTE])[kte]=(c1o1-q)/(c1o1+q)*(f_BW-f_TE+(f_BW+f_TE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BW+f_TE)-c6o1*c1o54*(-VeloX-VeloZ))/(c1o1+q) - c1o54 * drho;
//(D.f[dirTE])[kte]=zero;
f_TE_in=(c1o1-q)/(c1o1+q)*(f_BW-f_TE+(f_BW+f_TE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BW+f_TE))/(c1o1+q) - c1o54 * drho;
wallMomentumX -= f_BW+f_TE_in;
wallMomentumZ -= f_BW+f_TE_in;
}
q = q_dirBE[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o54* (drho/*+three*( vx1 -vx3)*/+c9o2*( vx1 -vx3)*( vx1 -vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirTW])[ktw]=(c1o1-q)/(c1o1+q)*(f_BE-f_TW+(f_BE+f_TW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BE+f_TW)-c6o1*c1o54*( VeloX-VeloZ))/(c1o1+q) - c1o54 * drho;
//(D.f[dirTW])[ktw]=zero;
f_TW_in=(c1o1-q)/(c1o1+q)*(f_BE-f_TW+(f_BE+f_TW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BE+f_TW))/(c1o1+q) - c1o54 * drho;
wallMomentumX += f_BE+f_TW_in;
wallMomentumZ -= f_BE+f_TW_in;
}
q = q_dirTW[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o54* (drho/*+three*(-vx1 +vx3)*/+c9o2*(-vx1 +vx3)*(-vx1 +vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirBE])[kbe]=(c1o1-q)/(c1o1+q)*(f_TW-f_BE+(f_TW+f_BE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TW+f_BE)-c6o1*c1o54*(-VeloX+VeloZ))/(c1o1+q) - c1o54 * drho;
//(D.f[dirBE])[kbe]=zero;
f_BE_in=(c1o1-q)/(c1o1+q)*(f_TW-f_BE+(f_TW+f_BE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TW+f_BE))/(c1o1+q) - c1o54 * drho;
wallMomentumX -= f_TW+f_BE_in;
wallMomentumZ += f_TW+f_BE_in;
}
q = q_dirTN[k];
if (q>=c0o1 && q<=c1o1)
{
{
feq=c1o54* (drho/*+three*( vx2+vx3)*/+c9o2*( vx2+vx3)*( vx2+vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirBS])[kbs]=(c1o1-q)/(c1o1+q)*(f_TN-f_BS+(f_TN+f_BS-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TN+f_BS)-c6o1*c1o54*( VeloY+VeloZ))/(c1o1+q) - c1o54 * drho;
//(D.f[dirBS])[kbs]=zero;
f_BS_in=(c1o1-q)/(c1o1+q)*(f_TN-f_BS+(f_TN+f_BS-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TN+f_BS))/(c1o1+q) - c1o54 * drho;
wallMomentumY += f_TN+f_BS_in;
wallMomentumZ += f_TN+f_BS_in;
}
q = q_dirBS[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o54* (drho/*+three*( -vx2-vx3)*/+c9o2*( -vx2-vx3)*( -vx2-vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirTN])[ktn]=(c1o1-q)/(c1o1+q)*(f_BS-f_TN+(f_BS+f_TN-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BS+f_TN)-c6o1*c1o54*( -VeloY-VeloZ))/(c1o1+q) - c1o54 * drho;
//(D.f[dirTN])[ktn]=zero;
f_TN_in=(c1o1-q)/(c1o1+q)*(f_BS-f_TN+(f_BS+f_TN-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BS+f_TN))/(c1o1+q) - c1o54 * drho;
wallMomentumY -= f_BS+f_TN_in;
wallMomentumZ -= f_BS+f_TN_in;
}
q = q_dirBN[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o54* (drho/*+three*( vx2-vx3)*/+c9o2*( vx2-vx3)*( vx2-vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirTS])[kts]=(c1o1-q)/(c1o1+q)*(f_BN-f_TS+(f_BN+f_TS-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BN+f_TS)-c6o1*c1o54*( VeloY-VeloZ))/(c1o1+q) - c1o54 * drho;
//(D.f[dirTS])[kts]=zero;
f_TS_in=(c1o1-q)/(c1o1+q)*(f_BN-f_TS+(f_BN+f_TS-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BN+f_TS))/(c1o1+q) - c1o54 * drho;
wallMomentumY += f_BN+f_TS_in;
wallMomentumZ -= f_BN+f_TS_in;
}
q = q_dirTS[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o54* (drho/*+three*( -vx2+vx3)*/+c9o2*( -vx2+vx3)*( -vx2+vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirBN])[kbn]=(c1o1-q)/(c1o1+q)*(f_TS-f_BN+(f_TS+f_BN-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TS+f_BN)-c6o1*c1o54*( -VeloY+VeloZ))/(c1o1+q) - c1o54 * drho;
//(D.f[dirBN])[kbn]=zero;
f_BN_in=(c1o1-q)/(c1o1+q)*(f_TS-f_BN+(f_TS+f_BN-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TS+f_BN))/(c1o1+q) - c1o54 * drho;
wallMomentumY -= f_TS+f_BN_in;
wallMomentumZ += f_TS+f_BN_in;
}
q = q_dirTNE[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o216*(drho/*+three*( vx1+vx2+vx3)*/+c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirBSW])[kbsw]=(c1o1-q)/(c1o1+q)*(f_TNE-f_BSW+(f_TNE+f_BSW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TNE+f_BSW)-c6o1*c1o216*( VeloX+VeloY+VeloZ))/(c1o1+q) - c1o216 * drho;
//(D.f[dirBSW])[kbsw]=zero;
f_BSW_in=(c1o1-q)/(c1o1+q)*(f_TNE-f_BSW+(f_TNE+f_BSW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TNE+f_BSW))/(c1o1+q) - c1o216 * drho;
wallMomentumX += f_TNE+f_BSW_in;
wallMomentumY += f_TNE+f_BSW_in;
wallMomentumZ += f_TNE+f_BSW_in;
}
q = q_dirBSW[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o216*(drho/*+three*(-vx1-vx2-vx3)*/+c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirTNE])[ktne]=(c1o1-q)/(c1o1+q)*(f_BSW-f_TNE+(f_BSW+f_TNE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BSW+f_TNE)-c6o1*c1o216*(-VeloX-VeloY-VeloZ))/(c1o1+q) - c1o216 * drho;
//(D.f[dirTNE])[ktne]=zero;
f_TNE_in=(c1o1-q)/(c1o1+q)*(f_BSW-f_TNE+(f_BSW+f_TNE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BSW+f_TNE))/(c1o1+q) - c1o216 * drho;
wallMomentumX -= f_BSW+f_TNE_in;
wallMomentumY -= f_BSW+f_TNE_in;
wallMomentumZ -= f_BSW+f_TNE_in;
}
q = q_dirBNE[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o216*(drho/*+three*( vx1+vx2-vx3)*/+c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirTSW])[ktsw]=(c1o1-q)/(c1o1+q)*(f_BNE-f_TSW+(f_BNE+f_TSW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BNE+f_TSW)-c6o1*c1o216*( VeloX+VeloY-VeloZ))/(c1o1+q) - c1o216 * drho;
//(D.f[dirTSW])[ktsw]=zero;
f_TSW_in=(c1o1-q)/(c1o1+q)*(f_BNE-f_TSW+(f_BNE+f_TSW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BNE+f_TSW))/(c1o1+q) - c1o216 * drho;
wallMomentumX += f_BNE+f_TSW_in;
wallMomentumY += f_BNE+f_TSW_in;
wallMomentumZ -= f_BNE+f_TSW_in;
}
q = q_dirTSW[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o216*(drho/*+three*(-vx1-vx2+vx3)*/+c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirBNE])[kbne]=(c1o1-q)/(c1o1+q)*(f_TSW-f_BNE+(f_TSW+f_BNE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TSW+f_BNE)-c6o1*c1o216*(-VeloX-VeloY+VeloZ))/(c1o1+q) - c1o216 * drho;
//(D.f[dirBNE])[kbne]=zero;
f_BNE_in=(c1o1-q)/(c1o1+q)*(f_TSW-f_BNE+(f_TSW+f_BNE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TSW+f_BNE))/(c1o1+q) - c1o216 * drho;
wallMomentumX -= f_TSW+f_BNE_in;
wallMomentumY -= f_TSW+f_BNE_in;
wallMomentumZ += f_TSW+f_BNE_in;
}
q = q_dirTSE[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o216*(drho/*+three*( vx1-vx2+vx3)*/+c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirBNW])[kbnw]=(c1o1-q)/(c1o1+q)*(f_TSE-f_BNW+(f_TSE+f_BNW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TSE+f_BNW)-c6o1*c1o216*( VeloX-VeloY+VeloZ))/(c1o1+q) - c1o216 * drho;
//(D.f[dirBNW])[kbnw]=zero;
f_BNW_in=(c1o1-q)/(c1o1+q)*(f_TSE-f_BNW+(f_TSE+f_BNW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TSE+f_BNW))/(c1o1+q) - c1o216 * drho;
wallMomentumX += f_TSE+f_BNW_in;
wallMomentumY -= f_TSE+f_BNW_in;
wallMomentumZ += f_TSE+f_BNW_in;
}
q = q_dirBNW[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o216*(drho/*+three*(-vx1+vx2-vx3)*/+c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirTSE])[ktse]=(c1o1-q)/(c1o1+q)*(f_BNW-f_TSE+(f_BNW+f_TSE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BNW+f_TSE)-c6o1*c1o216*(-VeloX+VeloY-VeloZ))/(c1o1+q) - c1o216 * drho;
//(D.f[dirTSE])[ktse]=zero;
f_TSE_in=(c1o1-q)/(c1o1+q)*(f_BNW-f_TSE+(f_BNW+f_TSE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BNW+f_TSE))/(c1o1+q) - c1o216 * drho;
wallMomentumX -= f_BNW+f_TSE_in;
wallMomentumY += f_BNW+f_TSE_in;
wallMomentumZ -= f_BNW+f_TSE_in;
}
q = q_dirBSE[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o216*(drho/*+three*( vx1-vx2-vx3)*/+c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirTNW])[ktnw]=(c1o1-q)/(c1o1+q)*(f_BSE-f_TNW+(f_BSE+f_TNW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BSE+f_TNW)-c6o1*c1o216*( VeloX-VeloY-VeloZ))/(c1o1+q) - c1o216 * drho;
//(D.f[dirTNW])[ktnw]=zero;
f_TNW_in=(c1o1-q)/(c1o1+q)*(f_BSE-f_TNW+(f_BSE+f_TNW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BSE+f_TNW))/(c1o1+q) - c1o216 * drho;
wallMomentumX += f_BSE+f_TNW_in;
wallMomentumY -= f_BSE+f_TNW_in;
wallMomentumZ -= f_BSE+f_TNW_in;
}
q = q_dirTNW[k];
if (q>=c0o1 && q<=c1o1)
{
feq=c1o216*(drho/*+three*(-vx1+vx2+vx3)*/+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3) * (c1o1 + drho)-cu_sq);
(D.f[dirBSE])[kbse]=(c1o1-q)/(c1o1+q)*(f_TNW-f_BSE+(f_TNW+f_BSE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TNW+f_BSE)-c6o1*c1o216*(-VeloX+VeloY+VeloZ))/(c1o1+q) - c1o216 * drho;
//(D.f[dirBSE])[kbse]=zero;
f_BSE_in=(c1o1-q)/(c1o1+q)*(f_TNW-f_BSE+(f_TNW+f_BSE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TNW+f_BSE))/(c1o1+q) - c1o216 * drho;
wallMomentumX -= f_TNW+f_BSE_in;
wallMomentumY += f_TNW+f_BSE_in;
wallMomentumZ += f_TNW+f_BSE_in;
}
// incoming f's from bounce back
// real f_E_in, f_W_in, f_N_in, f_S_in, f_T_in, f_B_in, f_NE_in, f_SW_in, f_SE_in, f_NW_in, f_TE_in, f_BW_in, f_BE_in,
// f_TW_in, f_TN_in, f_BS_in, f_BN_in, f_TS_in, f_TNE_in, f_TSW_in, f_TSE_in, f_TNW_in, f_BNE_in, f_BSW_in, f_BSE_in, f_BNW_in;
// // momentum exchanged with wall at rest
// real wallMomentumX = 0.0, wallMomentumY = 0.0, wallMomentumZ = 0.0;
// q = q_dirE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c2o27* (drho/*+three*( vx1 )*/+c9o2*( vx1 )*( vx1 ) * (c1o1 + drho)-cu_sq);
// f_W_in=(c1o1-q)/(c1o1+q)*(f_E-f_W+(f_E+f_W-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_E+f_W))/(c1o1+q) - c2o27 * drho;
// wallMomentumX += f_E+f_W_in;
// }
// q = q_dirW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c2o27* (drho/*+three*(-vx1 )*/+c9o2*(-vx1 )*(-vx1 ) * (c1o1 + drho)-cu_sq);
// f_E_in=(c1o1-q)/(c1o1+q)*(f_W-f_E+(f_W+f_E-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_W+f_E))/(c1o1+q) - c2o27 * drho;
// wallMomentumX -= f_W+f_E_in;
// }
// q = q_dirN[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c2o27* (drho/*+three*( vx2 )*/+c9o2*( vx2 )*( vx2 ) * (c1o1 + drho)-cu_sq);
// f_S_in=(c1o1-q)/(c1o1+q)*(f_N-f_S+(f_N+f_S-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_N+f_S))/(c1o1+q) - c2o27 * drho;
// wallMomentumY += f_N+f_S_in;
// }
// q = q_dirS[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c2o27* (drho/*+three*( -vx2 )*/+c9o2*( -vx2 )*( -vx2 ) * (c1o1 + drho)-cu_sq);
// f_N_in=(c1o1-q)/(c1o1+q)*(f_S-f_N+(f_S+f_N-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_S+f_N))/(c1o1+q) - c2o27 * drho;
// wallMomentumY -= f_S+f_N_in;
// }
// q = q_dirT[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c2o27* (drho/*+three*( vx3)*/+c9o2*( vx3)*( vx3) * (c1o1 + drho)-cu_sq);
// f_B_in=(c1o1-q)/(c1o1+q)*(f_T-f_B+(f_T+f_B-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_T+f_B))/(c1o1+q) - c2o27 * drho;
// wallMomentumZ += f_T+f_B_in;
// }
// q = q_dirB[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c2o27* (drho/*+three*( -vx3)*/+c9o2*( -vx3)*( -vx3) * (c1o1 + drho)-cu_sq);
// f_T_in=(c1o1-q)/(c1o1+q)*(f_B-f_T+(f_B+f_T-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_B+f_T))/(c1o1+q) - c2o27 * drho;
// wallMomentumZ -= f_B+f_T_in;
// }
// q = q_dirNE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o54* (drho/*+three*( vx1+vx2 )*/+c9o2*( vx1+vx2 )*( vx1+vx2 ) * (c1o1 + drho)-cu_sq);
// f_SW_in=(c1o1-q)/(c1o1+q)*(f_NE-f_SW+(f_NE+f_SW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_NE+f_SW))/(c1o1+q) - c1o54 * drho;
// wallMomentumX += f_NE+f_SW_in;
// wallMomentumY += f_NE+f_SW_in;
// }
// q = q_dirSW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o54* (drho/*+three*(-vx1-vx2 )*/+c9o2*(-vx1-vx2 )*(-vx1-vx2 ) * (c1o1 + drho)-cu_sq);
// f_NE_in=(c1o1-q)/(c1o1+q)*(f_SW-f_NE+(f_SW+f_NE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_SW+f_NE))/(c1o1+q) - c1o54 * drho;
// wallMomentumX -= f_SW+f_NE_in;
// wallMomentumY -= f_SW+f_NE_in;
// }
// q = q_dirSE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o54* (drho/*+three*( vx1-vx2 )*/+c9o2*( vx1-vx2 )*( vx1-vx2 ) * (c1o1 + drho)-cu_sq);
// f_NW_in=(c1o1-q)/(c1o1+q)*(f_SE-f_NW+(f_SE+f_NW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_SE+f_NW))/(c1o1+q) - c1o54 * drho;
// wallMomentumX += f_SE+f_NW_in;
// wallMomentumY -= f_SE+f_NW_in;
// }
// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// //Compute wall velocity
// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
real VeloX=0.0, VeloY=0.0, VeloZ=0.0;
// q = q_dirNW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o54* (drho/*+three*(-vx1+vx2 )*/+c9o2*(-vx1+vx2 )*(-vx1+vx2 ) * (c1o1 + drho)-cu_sq);
// f_SE_in=(c1o1-q)/(c1o1+q)*(f_NW-f_SE+(f_NW+f_SE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_NW+f_SE))/(c1o1+q) - c1o54 * drho;
// wallMomentumX -= f_NW+f_SE_in;
// wallMomentumY += f_NW+f_SE_in;
// }
// q = q_dirTE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o54* (drho/*+three*( vx1 +vx3)*/+c9o2*( vx1 +vx3)*( vx1 +vx3) * (c1o1 + drho)-cu_sq);
// f_BW_in=(c1o1-q)/(c1o1+q)*(f_TE-f_BW+(f_TE+f_BW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TE+f_BW))/(c1o1+q) - c1o54 * drho;
// wallMomentumX += f_TE+f_BW_in;
// wallMomentumZ += f_TE+f_BW_in;
// }
// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// //Add wall velocity and write f's
// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// q = q_dirBW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o54* (drho/*+three*(-vx1 -vx3)*/+c9o2*(-vx1 -vx3)*(-vx1 -vx3) * (c1o1 + drho)-cu_sq);
// f_TE_in=(c1o1-q)/(c1o1+q)*(f_BW-f_TE+(f_BW+f_TE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BW+f_TE))/(c1o1+q) - c1o54 * drho;
// wallMomentumX -= f_BW+f_TE_in;
// wallMomentumZ -= f_BW+f_TE_in;
// }
// q = q_dirBE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o54* (drho/*+three*( vx1 -vx3)*/+c9o2*( vx1 -vx3)*( vx1 -vx3) * (c1o1 + drho)-cu_sq);
// f_TW_in=(c1o1-q)/(c1o1+q)*(f_BE-f_TW+(f_BE+f_TW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BE+f_TW))/(c1o1+q) - c1o54 * drho;
// wallMomentumX += f_BE+f_TW_in;
// wallMomentumZ -= f_BE+f_TW_in;
// }
// q = q_dirTW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o54* (drho/*+three*(-vx1 +vx3)*/+c9o2*(-vx1 +vx3)*(-vx1 +vx3) * (c1o1 + drho)-cu_sq);
// f_BE_in=(c1o1-q)/(c1o1+q)*(f_TW-f_BE+(f_TW+f_BE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TW+f_BE))/(c1o1+q) - c1o54 * drho;
// wallMomentumX -= f_TW+f_BE_in;
// wallMomentumZ += f_TW+f_BE_in;
// }
// q = q_dirTN[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o54* (drho/*+three*( vx2+vx3)*/+c9o2*( vx2+vx3)*( vx2+vx3) * (c1o1 + drho)-cu_sq);
// f_BS_in=(c1o1-q)/(c1o1+q)*(f_TN-f_BS+(f_TN+f_BS-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TN+f_BS))/(c1o1+q) - c1o54 * drho;
// wallMomentumY += f_TN+f_BS_in;
// wallMomentumZ += f_TN+f_BS_in;
// }
// q = q_dirBS[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o54* (drho/*+three*( -vx2-vx3)*/+c9o2*( -vx2-vx3)*( -vx2-vx3) * (c1o1 + drho)-cu_sq);
// f_TN_in=(c1o1-q)/(c1o1+q)*(f_BS-f_TN+(f_BS+f_TN-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BS+f_TN))/(c1o1+q) - c1o54 * drho;
// wallMomentumY -= f_BS+f_TN_in;
// wallMomentumZ -= f_BS+f_TN_in;
// }
// q = q_dirBN[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o54* (drho/*+three*( vx2-vx3)*/+c9o2*( vx2-vx3)*( vx2-vx3) * (c1o1 + drho)-cu_sq);
// f_TS_in=(c1o1-q)/(c1o1+q)*(f_BN-f_TS+(f_BN+f_TS-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BN+f_TS))/(c1o1+q) - c1o54 * drho;
// wallMomentumY += f_BN+f_TS_in;
// wallMomentumZ -= f_BN+f_TS_in;
// }
// q = q_dirTS[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o54* (drho/*+three*( -vx2+vx3)*/+c9o2*( -vx2+vx3)*( -vx2+vx3) * (c1o1 + drho)-cu_sq);
// f_BN_in=(c1o1-q)/(c1o1+q)*(f_TS-f_BN+(f_TS+f_BN-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TS+f_BN))/(c1o1+q) - c1o54 * drho;
// wallMomentumY -= f_TS+f_BN_in;
// wallMomentumZ += f_TS+f_BN_in;
// }
// q = q_dirTNE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o216*(drho/*+three*( vx1+vx2+vx3)*/+c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3) * (c1o1 + drho)-cu_sq);
// f_BSW_in=(c1o1-q)/(c1o1+q)*(f_TNE-f_BSW+(f_TNE+f_BSW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TNE+f_BSW))/(c1o1+q) - c1o216 * drho;
// wallMomentumX += f_TNE+f_BSW_in;
// wallMomentumY += f_TNE+f_BSW_in;
// wallMomentumZ += f_TNE+f_BSW_in;
// }
// q = q_dirBSW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o216*(drho/*+three*(-vx1-vx2-vx3)*/+c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3) * (c1o1 + drho)-cu_sq);
// f_TNE_in=(c1o1-q)/(c1o1+q)*(f_BSW-f_TNE+(f_BSW+f_TNE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BSW+f_TNE))/(c1o1+q) - c1o216 * drho;
// wallMomentumX -= f_BSW+f_TNE_in;
// wallMomentumY -= f_BSW+f_TNE_in;
// wallMomentumZ -= f_BSW+f_TNE_in;
// }
// q = q_dirBNE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o216*(drho/*+three*( vx1+vx2-vx3)*/+c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3) * (c1o1 + drho)-cu_sq);
// f_TSW_in=(c1o1-q)/(c1o1+q)*(f_BNE-f_TSW+(f_BNE+f_TSW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BNE+f_TSW))/(c1o1+q) - c1o216 * drho;
// wallMomentumX += f_BNE+f_TSW_in;
// wallMomentumY += f_BNE+f_TSW_in;
// wallMomentumZ -= f_BNE+f_TSW_in;
// }
// q = q_dirTSW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o216*(drho/*+three*(-vx1-vx2+vx3)*/+c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3) * (c1o1 + drho)-cu_sq);
// f_BNE_in=(c1o1-q)/(c1o1+q)*(f_TSW-f_BNE+(f_TSW+f_BNE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TSW+f_BNE))/(c1o1+q) - c1o216 * drho;
// wallMomentumX -= f_TSW+f_BNE_in;
// wallMomentumY -= f_TSW+f_BNE_in;
// wallMomentumZ += f_TSW+f_BNE_in;
// }
// q = q_dirTSE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o216*(drho/*+three*( vx1-vx2+vx3)*/+c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3) * (c1o1 + drho)-cu_sq);
// f_BNW_in=(c1o1-q)/(c1o1+q)*(f_TSE-f_BNW+(f_TSE+f_BNW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TSE+f_BNW))/(c1o1+q) - c1o216 * drho;
// wallMomentumX += f_TSE+f_BNW_in;
// wallMomentumY -= f_TSE+f_BNW_in;
// wallMomentumZ += f_TSE+f_BNW_in;
// }
// q = q_dirBNW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o216*(drho/*+three*(-vx1+vx2-vx3)*/+c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3) * (c1o1 + drho)-cu_sq);
// f_TSE_in=(c1o1-q)/(c1o1+q)*(f_BNW-f_TSE+(f_BNW+f_TSE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BNW+f_TSE))/(c1o1+q) - c1o216 * drho;
// wallMomentumX -= f_BNW+f_TSE_in;
// wallMomentumY += f_BNW+f_TSE_in;
// wallMomentumZ -= f_BNW+f_TSE_in;
// }
// q = q_dirBSE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o216*(drho/*+three*( vx1-vx2-vx3)*/+c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3) * (c1o1 + drho)-cu_sq);
// f_TNW_in=(c1o1-q)/(c1o1+q)*(f_BSE-f_TNW+(f_BSE+f_TNW-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_BSE+f_TNW))/(c1o1+q) - c1o216 * drho;
// wallMomentumX += f_BSE+f_TNW_in;
// wallMomentumY -= f_BSE+f_TNW_in;
// wallMomentumZ -= f_BSE+f_TNW_in;
// }
// q = q_dirTNW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// feq=c1o216*(drho/*+three*(-vx1+vx2+vx3)*/+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3) * (c1o1 + drho)-cu_sq);
// f_BSE_in=(c1o1-q)/(c1o1+q)*(f_TNW-f_BSE+(f_TNW+f_BSE-c2o1*feq*om1)/(c1o1-om1))*c1o2+(q*(f_TNW+f_BSE))/(c1o1+q) - c1o216 * drho;
// wallMomentumX -= f_TNW+f_BSE_in;
// wallMomentumY += f_TNW+f_BSE_in;
// wallMomentumZ += f_TNW+f_BSE_in;
// }
// // ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// // //Compute wall velocity
// // ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// real VeloX=0.1, VeloY=0.0, VeloZ=0.0;
// // ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// // //Add wall velocity and write f's
// // ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// q = q_dirE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirW])[kw] = f_W_in - (c6o1*c2o27*( VeloX ))/(c1o1+q);
// }
q = q_dirE[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirW])[kw] = f_W_in - (c6o1*c2o27*( VeloX ))/(c1o1+q);
}
// q = q_dirW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirE])[ke] = f_E_in - (c6o1*c2o27*(-VeloX ))/(c1o1+q);
// }
q = q_dirW[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirE])[ke] = f_E_in - (c6o1*c2o27*(-VeloX ))/(c1o1+q);
}
// q = q_dirN[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirS])[ks] = f_S_in - (c6o1*c2o27*( VeloY ))/(c1o1+q);
// }
q = q_dirN[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirS])[ks] = f_S_in - (c6o1*c2o27*( VeloY ))/(c1o1+q);
}
// q = q_dirS[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirN])[kn] = f_N_in - (c6o1*c2o27*(-VeloY ))/(c1o1+q);
// }
q = q_dirS[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirN])[kn] = f_N_in - (c6o1*c2o27*(-VeloY ))/(c1o1+q);
}
// q = q_dirT[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirB])[kb] = f_B_in - (c6o1*c2o27*( VeloZ ))/(c1o1+q);
// }
q = q_dirT[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirB])[kb] = f_B_in - (c6o1*c2o27*( VeloZ ))/(c1o1+q);
}
// q = q_dirB[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirT])[kt] = f_T_in - (c6o1*c2o27*(-VeloZ ))/(c1o1+q);
// }
q = q_dirB[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirT])[kt] = f_T_in - (c6o1*c2o27*(-VeloZ ))/(c1o1+q);
}
// q = q_dirNE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirSW])[ksw] = f_SW_in - (c6o1*c1o54*(VeloX+VeloY))/(c1o1+q);
// }
q = q_dirNE[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirSW])[ksw] = f_SW_in - (c6o1*c1o54*(VeloX+VeloY))/(c1o1+q);
}
// q = q_dirSW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirNE])[kne] = f_NE_in - (c6o1*c1o54*(-VeloX-VeloY))/(c1o1+q);
// }
q = q_dirSW[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirNE])[kne] = f_NE_in - (c6o1*c1o54*(-VeloX-VeloY))/(c1o1+q);
}
// q = q_dirSE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirNW])[knw] = f_NW_in - (c6o1*c1o54*( VeloX-VeloY))/(c1o1+q);
// }
q = q_dirSE[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirNW])[knw] = f_NW_in - (c6o1*c1o54*( VeloX-VeloY))/(c1o1+q);
}
// q = q_dirNW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirSE])[kse] = f_SE_in - (c6o1*c1o54*(-VeloX+VeloY))/(c1o1+q);
// }
q = q_dirNW[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirSE])[kse] = f_SE_in - (c6o1*c1o54*(-VeloX+VeloY))/(c1o1+q);
}
// q = q_dirTE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirBW])[kbw] = f_BW_in - (c6o1*c1o54*( VeloX+VeloZ))/(c1o1+q);
// }
q = q_dirTE[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirBW])[kbw] = f_BW_in - (c6o1*c1o54*( VeloX+VeloZ))/(c1o1+q);
}
// q = q_dirBW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirTE])[kte] = f_TE_in - (c6o1*c1o54*(-VeloX-VeloZ))/(c1o1+q);
// }
q = q_dirBW[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirTE])[kte] = f_TE_in - (c6o1*c1o54*(-VeloX-VeloZ))/(c1o1+q);
}
// q = q_dirBE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirTW])[ktw] = f_TW_in - (c6o1*c1o54*( VeloX-VeloZ))/(c1o1+q);
// }
q = q_dirBE[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirTW])[ktw] = f_TW_in - (c6o1*c1o54*( VeloX-VeloZ))/(c1o1+q);
}
// q = q_dirTW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirBE])[kbe] = f_BE_in - (c6o1*c1o54*(-VeloX+VeloZ))/(c1o1+q);
// }
q = q_dirTW[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirBE])[kbe] = f_BE_in - (c6o1*c1o54*(-VeloX+VeloZ))/(c1o1+q);
}
// q = q_dirTN[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirBS])[kbs] = f_BS_in - (c6o1*c1o54*( VeloY+VeloZ))/(c1o1+q);
// }
q = q_dirTN[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirBS])[kbs] = f_BS_in - (c6o1*c1o54*( VeloY+VeloZ))/(c1o1+q);
}
// q = q_dirBS[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirTN])[ktn] = f_TN_in - (c6o1*c1o54*( -VeloY-VeloZ))/(c1o1+q);
// }
q = q_dirBS[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirTN])[ktn] = f_TN_in - (c6o1*c1o54*( -VeloY-VeloZ))/(c1o1+q);
}
// q = q_dirBN[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirTS])[kts] = f_TS_in - (c6o1*c1o54*( VeloY-VeloZ))/(c1o1+q);
// }
q = q_dirBN[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirTS])[kts] = f_TS_in - (c6o1*c1o54*( VeloY-VeloZ))/(c1o1+q);
}
// q = q_dirTS[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirBN])[kbn] = f_BN_in - (c6o1*c1o54*( -VeloY+VeloZ))/(c1o1+q);
// }
q = q_dirTS[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirBN])[kbn] = f_BN_in - (c6o1*c1o54*( -VeloY+VeloZ))/(c1o1+q);
}
// q = q_dirTNE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirBSW])[kbsw] = f_BSW_in - (c6o1*c1o216*( VeloX+VeloY+VeloZ))/(c1o1+q);
// }
q = q_dirTNE[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirBSW])[kbsw] = f_BSW_in - (c6o1*c1o216*( VeloX+VeloY+VeloZ))/(c1o1+q);
}
// q = q_dirBSW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirTNE])[ktne] = f_TNE_in - (c6o1*c1o216*(-VeloX-VeloY-VeloZ))/(c1o1+q);
// }
q = q_dirBSW[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirTNE])[ktne] = f_TNE_in - (c6o1*c1o216*(-VeloX-VeloY-VeloZ))/(c1o1+q);
}
// q = q_dirBNE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirTSW])[ktsw] = f_TSW_in - (c6o1*c1o216*( VeloX+VeloY-VeloZ))/(c1o1+q);
// }
q = q_dirBNE[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirTSW])[ktsw] = f_TSW_in - (c6o1*c1o216*( VeloX+VeloY-VeloZ))/(c1o1+q);
}
// q = q_dirTSW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirBNE])[kbne] = f_BNE_in - (c6o1*c1o216*(-VeloX-VeloY+VeloZ))/(c1o1+q);
// }
q = q_dirTSW[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirBNE])[kbne] = f_BNE_in - (c6o1*c1o216*(-VeloX-VeloY+VeloZ))/(c1o1+q);
}
// q = q_dirTSE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirBNW])[kbnw] = f_BNW_in - (c6o1*c1o216*( VeloX-VeloY+VeloZ))/(c1o1+q);
// }
q = q_dirTSE[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirBNW])[kbnw] = f_BNW_in - (c6o1*c1o216*( VeloX-VeloY+VeloZ))/(c1o1+q);
}
// q = q_dirBNW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirTSE])[ktse] = f_TSE_in - (c6o1*c1o216*(-VeloX+VeloY-VeloZ))/(c1o1+q);
// }
q = q_dirBNW[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirTSE])[ktse] = f_TSE_in - (c6o1*c1o216*(-VeloX+VeloY-VeloZ))/(c1o1+q);
}
// q = q_dirBSE[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirTNW])[ktnw] = f_TNW_in - (c6o1*c1o216*( VeloX-VeloY-VeloZ))/(c1o1+q);
// }
q = q_dirBSE[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirTNW])[ktnw] = f_TNW_in - (c6o1*c1o216*( VeloX-VeloY-VeloZ))/(c1o1+q);
}
// q = q_dirTNW[k];
// if (q>=c0o1 && q<=c1o1)
// {
// (D.f[dirBSE])[kbse] = f_BSE_in - (c6o1*c1o216*(-VeloX+VeloY+VeloZ))/(c1o1+q);
// }
q = q_dirTNW[k];
if (q>=c0o1 && q<=c1o1)
{
(D.f[dirBSE])[kbse] = f_BSE_in - (c6o1*c1o216*(-VeloX+VeloY+VeloZ))/(c1o1+q);
}
}
}
\ No newline at end of file
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment