diff --git a/src/gpu/VirtualFluids_GPU/AdvectionDiffusion/AdvectionDiffusion.cpp b/src/gpu/VirtualFluids_GPU/AdvectionDiffusion/AdvectionDiffusion.cpp index 3c408bd657fe33c48a565d6a0b8ec39bdad978a6..f466cb5c0a21fdf7119cf7bc682d2b432a5e3a96 100644 --- a/src/gpu/VirtualFluids_GPU/AdvectionDiffusion/AdvectionDiffusion.cpp +++ b/src/gpu/VirtualFluids_GPU/AdvectionDiffusion/AdvectionDiffusion.cpp @@ -125,19 +125,6 @@ void calcAD(SPtr<Parameter> para) para->getParD()->forcing, para->getParD()->isEvenTimestep); - /*CentralMomentsAdvectionDiffusionDeviceKernel( - para->getParD()->numberofthreads, - para->getParD()->omegaD, - para->getParD()->typeOfGridNode, - para->getParD()->neighborX, - para->getParD()->neighborY, - para->getParD()->neighborZ, - para->getParD()->distributions.f[0], - para->getParD()->distributionsAD.f[0], - para->getParD()->numberOfNodes, - para->getParD()->forcing, - para->getParD()->isEvenTimestep);*/ - if (para->getParD()->numberOfSlipBCnodes > 1) { ADSlipVelDevComp( para->getParD()->numberofthreads, diff --git a/src/gpu/VirtualFluids_GPU/GPU/AdvectionDiffusion27chim.cu b/src/gpu/VirtualFluids_GPU/GPU/AdvectionDiffusion27chim.cu index 830705c89426ae19216e73252ee8d6e9e0e2b5bc..4a593c06dd7a71f33b29af9d138eec892a4538d2 100644 --- a/src/gpu/VirtualFluids_GPU/GPU/AdvectionDiffusion27chim.cu +++ b/src/gpu/VirtualFluids_GPU/GPU/AdvectionDiffusion27chim.cu @@ -35,36 +35,6 @@ #include "LBM/D3Q27.h" #include "Core/RealConstants.h" -//////////////////////////////////////////////////////////////////////////////// -//! \brief forward chimera transformation \ref forwardInverseChimeraWithK -//! Transformation from distributions to central moments according to Eq. (6)-(14) in \ref -//! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a> -//! Modified for lower round-off errors. -inline __device__ void forwardInverseChimeraWithKincompressible(real &mfa, real &mfb, real &mfc, real vv, real v2, real Kinverse, real K, real oneMinusRho) { - real m2 = mfa + mfc; - real m1 = mfc - mfa; - real m0 = m2 + mfb; - mfa = m0; - m0 *= Kinverse; - m0 += oneMinusRho; - mfb = (m1*Kinverse - m0 * vv) * K; - mfc = ((m2 - c2o1* m1 * vv)*Kinverse + v2 * m0) * K; -} - -//////////////////////////////////////////////////////////////////////////////// -//! \brief backward chimera transformation \ref backwardInverseChimeraWithK -//! Transformation from central moments to distributions according to Eq. (57)-(65) in \ref -//! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a> -//! Modified for lower round-off errors. -inline __device__ void backwardInverseChimeraWithKincompressible(real &mfa, real &mfb, real &mfc, real vv, real v2, real Kinverse, real K, real oneMinusRho) { - real m0 = (((mfc - mfb) * c1o2 + mfb * vv)*Kinverse + (mfa*Kinverse + oneMinusRho) * (v2 - vv) * c1o2) * K; - real m1 = (((mfa - mfc) - c2o1 * mfb * vv)*Kinverse + (mfa*Kinverse + oneMinusRho) * (-v2)) * K; - mfc = (((mfc + mfb) * c1o2 + mfb * vv)*Kinverse + (mfa*Kinverse + oneMinusRho) * (v2 + vv) * c1o2) * K; - mfa = m0; - mfb = m1; -} - - //////////////////////////////////////////////////////////////////////////////// //! \brief forward chimera transformation \ref forwardChimera @@ -561,514 +531,3 @@ extern "C" __global__ void Factorized_Central_Moments_Advection_Diffusion_Device } //////////////////////////////////////////////////////////////////////////////// - - - - - - - - - - - - - - -//////////////////////////////////////////////////////////////////////////////// -extern "C" __global__ void Central_Moments_Advection_Diffusion_Device_Kernel( - real omegaD, - uint* typeOfGridNode, - uint* neighborX, - uint* neighborY, - uint* neighborZ, - real* distributions, - real* distributionsAD, - int size_Mat, - real* forces, - bool isEvenTimestep) -{ - ////////////////////////////////////////////////////////////////////////// - //! Advection Diffusion Kernel based on Central Moments - //! - //! to be executed before Fluid Kernel - //! -> wrong advection direction if executed in wrong order - //! - //////////////////////////////////////////////////////////////////////////////// - //! - Get node index coordinates from thredIdx, blockIdx, blockDim and gridDim. - //! - const unsigned x = threadIdx.x; - const unsigned y = blockIdx.x; - const unsigned z = blockIdx.y; - - const unsigned nx = blockDim.x; - const unsigned ny = gridDim.x; - - const unsigned k = nx*(ny*z + y) + x; - - ////////////////////////////////////////////////////////////////////////// - // run for all indices in size_Mat and fluid nodes - if ((k < size_Mat) && (typeOfGridNode[k] == GEO_FLUID)) - { - ////////////////////////////////////////////////////////////////////////// - //! - Read distributions: style of reading and writing the distributions from/to stored arrays dependent on timestep is based on the esoteric twist algorithm \ref - //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a> - //! - Distributions27 dist; - if (isEvenTimestep) - { - dist.f[dirE] = &distributions[dirE *size_Mat]; - dist.f[dirW] = &distributions[dirW *size_Mat]; - dist.f[dirN] = &distributions[dirN *size_Mat]; - dist.f[dirS] = &distributions[dirS *size_Mat]; - dist.f[dirT] = &distributions[dirT *size_Mat]; - dist.f[dirB] = &distributions[dirB *size_Mat]; - dist.f[dirNE] = &distributions[dirNE *size_Mat]; - dist.f[dirSW] = &distributions[dirSW *size_Mat]; - dist.f[dirSE] = &distributions[dirSE *size_Mat]; - dist.f[dirNW] = &distributions[dirNW *size_Mat]; - dist.f[dirTE] = &distributions[dirTE *size_Mat]; - dist.f[dirBW] = &distributions[dirBW *size_Mat]; - dist.f[dirBE] = &distributions[dirBE *size_Mat]; - dist.f[dirTW] = &distributions[dirTW *size_Mat]; - dist.f[dirTN] = &distributions[dirTN *size_Mat]; - dist.f[dirBS] = &distributions[dirBS *size_Mat]; - dist.f[dirBN] = &distributions[dirBN *size_Mat]; - dist.f[dirTS] = &distributions[dirTS *size_Mat]; - dist.f[dirREST] = &distributions[dirREST*size_Mat]; - dist.f[dirTNE] = &distributions[dirTNE *size_Mat]; - dist.f[dirTSW] = &distributions[dirTSW *size_Mat]; - dist.f[dirTSE] = &distributions[dirTSE *size_Mat]; - dist.f[dirTNW] = &distributions[dirTNW *size_Mat]; - dist.f[dirBNE] = &distributions[dirBNE *size_Mat]; - dist.f[dirBSW] = &distributions[dirBSW *size_Mat]; - dist.f[dirBSE] = &distributions[dirBSE *size_Mat]; - dist.f[dirBNW] = &distributions[dirBNW *size_Mat]; - } - else - { - dist.f[dirW] = &distributions[dirE *size_Mat]; - dist.f[dirE] = &distributions[dirW *size_Mat]; - dist.f[dirS] = &distributions[dirN *size_Mat]; - dist.f[dirN] = &distributions[dirS *size_Mat]; - dist.f[dirB] = &distributions[dirT *size_Mat]; - dist.f[dirT] = &distributions[dirB *size_Mat]; - dist.f[dirSW] = &distributions[dirNE *size_Mat]; - dist.f[dirNE] = &distributions[dirSW *size_Mat]; - dist.f[dirNW] = &distributions[dirSE *size_Mat]; - dist.f[dirSE] = &distributions[dirNW *size_Mat]; - dist.f[dirBW] = &distributions[dirTE *size_Mat]; - dist.f[dirTE] = &distributions[dirBW *size_Mat]; - dist.f[dirTW] = &distributions[dirBE *size_Mat]; - dist.f[dirBE] = &distributions[dirTW *size_Mat]; - dist.f[dirBS] = &distributions[dirTN *size_Mat]; - dist.f[dirTN] = &distributions[dirBS *size_Mat]; - dist.f[dirTS] = &distributions[dirBN *size_Mat]; - dist.f[dirBN] = &distributions[dirTS *size_Mat]; - dist.f[dirREST] = &distributions[dirREST*size_Mat]; - dist.f[dirBSW] = &distributions[dirTNE *size_Mat]; - dist.f[dirBNE] = &distributions[dirTSW *size_Mat]; - dist.f[dirBNW] = &distributions[dirTSE *size_Mat]; - dist.f[dirBSE] = &distributions[dirTNW *size_Mat]; - dist.f[dirTSW] = &distributions[dirBNE *size_Mat]; - dist.f[dirTNE] = &distributions[dirBSW *size_Mat]; - dist.f[dirTNW] = &distributions[dirBSE *size_Mat]; - dist.f[dirTSE] = &distributions[dirBNW *size_Mat]; - } - //////////////////////////////////////////////////////////////////////////////// - Distributions27 distAD; - if (isEvenTimestep) - { - distAD.f[dirE] = &distributionsAD[dirE *size_Mat]; - distAD.f[dirW] = &distributionsAD[dirW *size_Mat]; - distAD.f[dirN] = &distributionsAD[dirN *size_Mat]; - distAD.f[dirS] = &distributionsAD[dirS *size_Mat]; - distAD.f[dirT] = &distributionsAD[dirT *size_Mat]; - distAD.f[dirB] = &distributionsAD[dirB *size_Mat]; - distAD.f[dirNE] = &distributionsAD[dirNE *size_Mat]; - distAD.f[dirSW] = &distributionsAD[dirSW *size_Mat]; - distAD.f[dirSE] = &distributionsAD[dirSE *size_Mat]; - distAD.f[dirNW] = &distributionsAD[dirNW *size_Mat]; - distAD.f[dirTE] = &distributionsAD[dirTE *size_Mat]; - distAD.f[dirBW] = &distributionsAD[dirBW *size_Mat]; - distAD.f[dirBE] = &distributionsAD[dirBE *size_Mat]; - distAD.f[dirTW] = &distributionsAD[dirTW *size_Mat]; - distAD.f[dirTN] = &distributionsAD[dirTN *size_Mat]; - distAD.f[dirBS] = &distributionsAD[dirBS *size_Mat]; - distAD.f[dirBN] = &distributionsAD[dirBN *size_Mat]; - distAD.f[dirTS] = &distributionsAD[dirTS *size_Mat]; - distAD.f[dirREST] = &distributionsAD[dirREST*size_Mat]; - distAD.f[dirTNE] = &distributionsAD[dirTNE *size_Mat]; - distAD.f[dirTSW] = &distributionsAD[dirTSW *size_Mat]; - distAD.f[dirTSE] = &distributionsAD[dirTSE *size_Mat]; - distAD.f[dirTNW] = &distributionsAD[dirTNW *size_Mat]; - distAD.f[dirBNE] = &distributionsAD[dirBNE *size_Mat]; - distAD.f[dirBSW] = &distributionsAD[dirBSW *size_Mat]; - distAD.f[dirBSE] = &distributionsAD[dirBSE *size_Mat]; - distAD.f[dirBNW] = &distributionsAD[dirBNW *size_Mat]; - } - else - { - distAD.f[dirW] = &distributionsAD[dirE *size_Mat]; - distAD.f[dirE] = &distributionsAD[dirW *size_Mat]; - distAD.f[dirS] = &distributionsAD[dirN *size_Mat]; - distAD.f[dirN] = &distributionsAD[dirS *size_Mat]; - distAD.f[dirB] = &distributionsAD[dirT *size_Mat]; - distAD.f[dirT] = &distributionsAD[dirB *size_Mat]; - distAD.f[dirSW] = &distributionsAD[dirNE *size_Mat]; - distAD.f[dirNE] = &distributionsAD[dirSW *size_Mat]; - distAD.f[dirNW] = &distributionsAD[dirSE *size_Mat]; - distAD.f[dirSE] = &distributionsAD[dirNW *size_Mat]; - distAD.f[dirBW] = &distributionsAD[dirTE *size_Mat]; - distAD.f[dirTE] = &distributionsAD[dirBW *size_Mat]; - distAD.f[dirTW] = &distributionsAD[dirBE *size_Mat]; - distAD.f[dirBE] = &distributionsAD[dirTW *size_Mat]; - distAD.f[dirBS] = &distributionsAD[dirTN *size_Mat]; - distAD.f[dirTN] = &distributionsAD[dirBS *size_Mat]; - distAD.f[dirTS] = &distributionsAD[dirBN *size_Mat]; - distAD.f[dirBN] = &distributionsAD[dirTS *size_Mat]; - distAD.f[dirREST] = &distributionsAD[dirREST*size_Mat]; - distAD.f[dirBSW] = &distributionsAD[dirTNE *size_Mat]; - distAD.f[dirBNE] = &distributionsAD[dirTSW *size_Mat]; - distAD.f[dirBNW] = &distributionsAD[dirTSE *size_Mat]; - distAD.f[dirBSE] = &distributionsAD[dirTNW *size_Mat]; - distAD.f[dirTSW] = &distributionsAD[dirBNE *size_Mat]; - distAD.f[dirTNE] = &distributionsAD[dirBSW *size_Mat]; - distAD.f[dirTNW] = &distributionsAD[dirBSE *size_Mat]; - distAD.f[dirTSE] = &distributionsAD[dirBNW *size_Mat]; - } - //////////////////////////////////////////////////////////////////////////////// - //! - Set neighbor indices (necessary for indirect addressing) - uint kw = neighborX[k]; - uint ks = neighborY[k]; - uint kb = neighborZ[k]; - uint ksw = neighborY[kw]; - uint kbw = neighborZ[kw]; - uint kbs = neighborZ[ks]; - uint kbsw = neighborZ[ksw]; - //////////////////////////////////////////////////////////////////////////////////// - //! - Set local distributions Fluid - //! - real fcbb = (dist.f[dirE])[k]; - real fabb = (dist.f[dirW])[kw]; - real fbcb = (dist.f[dirN])[k]; - real fbab = (dist.f[dirS])[ks]; - real fbbc = (dist.f[dirT])[k]; - real fbba = (dist.f[dirB])[kb]; - real fccb = (dist.f[dirNE])[k]; - real faab = (dist.f[dirSW])[ksw]; - real fcab = (dist.f[dirSE])[ks]; - real facb = (dist.f[dirNW])[kw]; - real fcbc = (dist.f[dirTE])[k]; - real faba = (dist.f[dirBW])[kbw]; - real fcba = (dist.f[dirBE])[kb]; - real fabc = (dist.f[dirTW])[kw]; - real fbcc = (dist.f[dirTN])[k]; - real fbaa = (dist.f[dirBS])[kbs]; - real fbca = (dist.f[dirBN])[kb]; - real fbac = (dist.f[dirTS])[ks]; - real fbbb = (dist.f[dirREST])[k]; - real fccc = (dist.f[dirTNE])[k]; - real faac = (dist.f[dirTSW])[ksw]; - real fcac = (dist.f[dirTSE])[ks]; - real facc = (dist.f[dirTNW])[kw]; - real fcca = (dist.f[dirBNE])[kb]; - real faaa = (dist.f[dirBSW])[kbsw]; - real fcaa = (dist.f[dirBSE])[kbs]; - real faca = (dist.f[dirBNW])[kbw]; - //////////////////////////////////////////////////////////////////////////////////// - //! - Set local distributions Advection Diffusion - //! - real mfcbb = (distAD.f[dirE])[k]; - real mfabb = (distAD.f[dirW])[kw]; - real mfbcb = (distAD.f[dirN])[k]; - real mfbab = (distAD.f[dirS])[ks]; - real mfbbc = (distAD.f[dirT])[k]; - real mfbba = (distAD.f[dirB])[kb]; - real mfccb = (distAD.f[dirNE])[k]; - real mfaab = (distAD.f[dirSW])[ksw]; - real mfcab = (distAD.f[dirSE])[ks]; - real mfacb = (distAD.f[dirNW])[kw]; - real mfcbc = (distAD.f[dirTE])[k]; - real mfaba = (distAD.f[dirBW])[kbw]; - real mfcba = (distAD.f[dirBE])[kb]; - real mfabc = (distAD.f[dirTW])[kw]; - real mfbcc = (distAD.f[dirTN])[k]; - real mfbaa = (distAD.f[dirBS])[kbs]; - real mfbca = (distAD.f[dirBN])[kb]; - real mfbac = (distAD.f[dirTS])[ks]; - real mfbbb = (distAD.f[dirREST])[k]; - real mfccc = (distAD.f[dirTNE])[k]; - real mfaac = (distAD.f[dirTSW])[ksw]; - real mfcac = (distAD.f[dirTSE])[ks]; - real mfacc = (distAD.f[dirTNW])[kw]; - real mfcca = (distAD.f[dirBNE])[kb]; - real mfaaa = (distAD.f[dirBSW])[kbsw]; - real mfcaa = (distAD.f[dirBSE])[kbs]; - real mfaca = (distAD.f[dirBNW])[kbw]; - //////////////////////////////////////////////////////////////////////////////////// - //! - Calculate density and velocity using pyramid summation for low round-off errors as in Eq. (J1)-(J3) \ref - //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa.2015.05.001 ]</b></a> - //! - //////////////////////////////////////////////////////////////////////////////////// - // fluid component - real drhoFluid = - ((((fccc + faaa) + (faca + fcac)) + ((facc + fcaa) + (faac + fcca))) + - (((fbac + fbca) + (fbaa + fbcc)) + ((fabc + fcba) + (faba + fcbc)) + ((facb + fcab) + (faab + fccb))) + - ((fabb + fcbb) + (fbab + fbcb) + (fbba + fbbc))) + fbbb; - - real rhoFluid = c1o1 + drhoFluid; - real OOrhoFluid = c1o1 / rhoFluid; - - real vvx = - ((((fccc - faaa) + (fcac - faca)) + ((fcaa - facc) + (fcca - faac))) + - (((fcba - fabc) + (fcbc - faba)) + ((fcab - facb) + (fccb - faab))) + - (fcbb - fabb)) * OOrhoFluid; - real vvy = - ((((fccc - faaa) + (faca - fcac)) + ((facc - fcaa) + (fcca - faac))) + - (((fbca - fbac) + (fbcc - fbaa)) + ((facb - fcab) + (fccb - faab))) + - (fbcb - fbab)) * OOrhoFluid; - real vvz = - ((((fccc - faaa) + (fcac - faca)) + ((facc - fcaa) + (faac - fcca))) + - (((fbac - fbca) + (fbcc - fbaa)) + ((fabc - fcba) + (fcbc - faba))) + - (fbbc - fbba)) * OOrhoFluid; - //////////////////////////////////////////////////////////////////////////////////// - // second component - real rho = - ((((mfccc + mfaaa) + (mfaca + mfcac)) + ((mfacc + mfcaa) + (mfaac + mfcca))) + - (((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) + ((mfacb + mfcab) + (mfaab + mfccb))) + - ((mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc))) + mfbbb; - //////////////////////////////////////////////////////////////////////////////////// - //! - Add half of the acceleration (body force) to the velocity as in Eq. (42) \ref - //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa.2015.05.001 ]</b></a> - //! - real fx = forces[0]; - real fy = forces[1]; - real fz = -rho*forces[2]; - vvx += fx*c1o2; - vvy += fy*c1o2; - vvz += fz*c1o2; - //////////////////////////////////////////////////////////////////////////////////// - real oneMinusRho = c1o1 - rho; - - real cx = - ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) + - (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) + - (mfcbb - mfabb)); - real cy = - ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) + - (((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) + - (mfbcb - mfbab)); - real cz = - ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) + - (((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) + - (mfbbc - mfbba)); - - //////////////////////////////////////////////////////////////////////////////////// - // calculate the square of velocities for this lattice node - real cx2 = cx*cx; - real cy2 = cy*cy; - real cz2 = cz*cz; - //////////////////////////////////////////////////////////////////////////////////// - //! - Chimera transform from well conditioned distributions to central moments as defined in Appendix J in \ref - //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa.2015.05.001 ]</b></a> - //! see also Eq. (6)-(14) in \ref - //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a> - //! - //////////////////////////////////////////////////////////////////////////////////// - // Z - Dir - forwardInverseChimeraWithKincompressible(mfaaa, mfaab, mfaac, cz, cz2, c36o1, c1o36, oneMinusRho); - forwardInverseChimeraWithKincompressible(mfaba, mfabb, mfabc, cz, cz2, c9o1, c1o9 , oneMinusRho); - forwardInverseChimeraWithKincompressible(mfaca, mfacb, mfacc, cz, cz2, c36o1, c1o36, oneMinusRho); - forwardInverseChimeraWithKincompressible(mfbaa, mfbab, mfbac, cz, cz2, c9o1, c1o9 , oneMinusRho); - forwardInverseChimeraWithKincompressible(mfbba, mfbbb, mfbbc, cz, cz2, c9o4, c4o9 , oneMinusRho); - forwardInverseChimeraWithKincompressible(mfbca, mfbcb, mfbcc, cz, cz2, c9o1, c1o9 , oneMinusRho); - forwardInverseChimeraWithKincompressible(mfcaa, mfcab, mfcac, cz, cz2, c36o1, c1o36, oneMinusRho); - forwardInverseChimeraWithKincompressible(mfcba, mfcbb, mfcbc, cz, cz2, c9o1, c1o9 , oneMinusRho); - forwardInverseChimeraWithKincompressible(mfcca, mfccb, mfccc, cz, cz2, c36o1, c1o36, oneMinusRho); - - //////////////////////////////////////////////////////////////////////////////////// - // Y - Dir - forwardInverseChimeraWithKincompressible(mfaaa, mfaba, mfaca, cy, cy2, c6o1, c1o6 , oneMinusRho); - forwardChimera( mfaab, mfabb, mfacb, cy, cy2); - forwardInverseChimeraWithKincompressible(mfaac, mfabc, mfacc, cy, cy2, c18o1, c1o18, oneMinusRho); - forwardInverseChimeraWithKincompressible(mfbaa, mfbba, mfbca, cy, cy2, c3o2, c2o3 , oneMinusRho); - forwardChimera( mfbab, mfbbb, mfbcb, cy, cy2); - forwardInverseChimeraWithKincompressible(mfbac, mfbbc, mfbcc, cy, cy2, c9o2, c2o9 , oneMinusRho); - forwardInverseChimeraWithKincompressible(mfcaa, mfcba, mfcca, cy, cy2, c6o1, c1o6 , oneMinusRho); - forwardChimera( mfcab, mfcbb, mfccb, cy, cy2); - forwardInverseChimeraWithKincompressible(mfcac, mfcbc, mfccc, cy, cy2, c18o1, c1o18, oneMinusRho); - - //////////////////////////////////////////////////////////////////////////////////// - // X - Dir - forwardInverseChimeraWithKincompressible(mfaaa, mfbaa, mfcaa, cx, cx2, c1o1, c1o1, oneMinusRho); - forwardChimera( mfaba, mfbba, mfcba, cx, cx2); - forwardInverseChimeraWithKincompressible(mfaca, mfbca, mfcca, cx, cx2, c3o1, c1o3, oneMinusRho); - forwardChimera( mfaab, mfbab, mfcab, cx, cx2); - forwardChimera( mfabb, mfbbb, mfcbb, cx, cx2); - forwardChimera( mfacb, mfbcb, mfccb, cx, cx2); - forwardInverseChimeraWithKincompressible(mfaac, mfbac, mfcac, cx, cx2, c3o1, c1o3, oneMinusRho); - forwardChimera( mfabc, mfbbc, mfcbc, cx, cx2); - forwardInverseChimeraWithKincompressible(mfacc, mfbcc, mfccc, cx, cx2, c3o1, c1o9, oneMinusRho); - - //////////////////////////////////////////////////////////////////////////////////// - //! - Factorized central moments for Advection Diffusion Equation - Eq. (15)-(16) in \ref - //! <a href="https://doi.org/10.1016/j.advwatres.2015.09.015"><b>[ X. Yang et al. (2016), DOI: 10.1016/j.advwatres.2015.09.015]</b></a> - //! - - // linearized orthogonalization of 3rd order central moments - real Mabc = mfabc - mfaba*c1o3; - real Mbca = mfbca - mfbaa*c1o3; - real Macb = mfacb - mfaab*c1o3; - real Mcba = mfcba - mfaba*c1o3; - real Mcab = mfcab - mfaab*c1o3; - real Mbac = mfbac - mfbaa*c1o3; - // linearized orthogonalization of 5th order central moments - real Mcbc = mfcbc - mfaba*c1o9; - real Mbcc = mfbcc - mfbaa*c1o9; - real Mccb = mfccb - mfaab*c1o9; - - // collision of 1st order moments - cx = cx * (c1o1 - omegaD) + omegaD * vvx * rho; - cy = cy * (c1o1 - omegaD) + omegaD * vvy * rho; - cz = cz * (c1o1 - omegaD) + omegaD * vvz * rho; - - cx2 = cx*cx; - cy2 = cy*cy; - cz2 = cz*cz; - - // equilibration of 2nd order moments - mfbba = c0o1; - mfbab = c0o1; - mfabb = c0o1; - - mfcaa = c1o3 * rho; - mfaca = c1o3 * rho; - mfaac = c1o3 * rho; - - - //real omega2 = 1.0f;// omegaD; - //mfbba *= (c1o1 - omega2); - //mfbab *= (c1o1 - omega2); - //mfabb *= (c1o1 - omega2); - - //mfcaa = mfcaa*(c1o1 - omega2) + omega2*c1o3 * rho; - //mfaca = mfaca*(c1o1 - omega2) + omega2*c1o3 * rho; - //mfaac = mfaac*(c1o1 - omega2) + omega2*c1o3 * rho; - - // equilibration of 3rd order moments - Mabc = c0o1; - Mbca = c0o1; - Macb = c0o1; - Mcba = c0o1; - Mcab = c0o1; - Mbac = c0o1; - mfbbb = c0o1; - - // from linearized orthogonalization 3rd order central moments to central moments - mfabc = Mabc + mfaba*c1o3; - mfbca = Mbca + mfbaa*c1o3; - mfacb = Macb + mfaab*c1o3; - mfcba = Mcba + mfaba*c1o3; - mfcab = Mcab + mfaab*c1o3; - mfbac = Mbac + mfbaa*c1o3; - - // equilibration of 4th order moments - mfacc = c1o9 * rho; - mfcac = c1o9 * rho; - mfcca = c1o9 * rho; - - mfcbb = c0o1; - mfbcb = c0o1; - mfbbc = c0o1; - - // equilibration of 5th order moments - Mcbc = c0o1; - Mbcc = c0o1; - Mccb = c0o1; - - // from linearized orthogonalization 5th order central moments to central moments - mfcbc = Mcbc + mfaba*c1o9; - mfbcc = Mbcc + mfbaa*c1o9; - mfccb = Mccb + mfaab*c1o9; - - // equilibration of 6th order moment - mfccc = c1o27 * rho; - - //////////////////////////////////////////////////////////////////////////////////// - //! - Chimera transform from central moments to well conditioned distributions as defined in Appendix J in - //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa.2015.05.001 ]</b></a> - //! see also Eq. (88)-(96) in - //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a> - //! - //////////////////////////////////////////////////////////////////////////////////// - // X - Dir - backwardInverseChimeraWithKincompressible(mfaaa, mfbaa, mfcaa, cx, cx2, c1o1, c1o1 , oneMinusRho); - backwardChimera( mfaba, mfbba, mfcba, cx, cx2); - backwardInverseChimeraWithKincompressible(mfaca, mfbca, mfcca, cx, cx2, c3o1, c1o3 , oneMinusRho); - backwardChimera( mfaab, mfbab, mfcab, cx, cx2); - backwardChimera( mfabb, mfbbb, mfcbb, cx, cx2); - backwardChimera( mfacb, mfbcb, mfccb, cx, cx2); - backwardInverseChimeraWithKincompressible(mfaac, mfbac, mfcac, cx, cx2, c3o1, c1o3 , oneMinusRho); - backwardChimera( mfabc, mfbbc, mfcbc, cx, cx2); - backwardInverseChimeraWithKincompressible(mfacc, mfbcc, mfccc, cx, cx2, c9o1, c1o9 , oneMinusRho); - - //////////////////////////////////////////////////////////////////////////////////// - // Y - Dir - backwardInverseChimeraWithKincompressible(mfaaa, mfaba, mfaca, cy, cy2, c6o1 , c1o6 , oneMinusRho); - backwardChimera( mfaab, mfabb, mfacb, cy, cy2); - backwardInverseChimeraWithKincompressible(mfaac, mfabc, mfacc, cy, cy2, c18o1, c1o18, oneMinusRho); - backwardInverseChimeraWithKincompressible(mfbaa, mfbba, mfbca, cy, cy2, c3o2 , c2o3 , oneMinusRho); - backwardChimera( mfbab, mfbbb, mfbcb, cy, cy2); - backwardInverseChimeraWithKincompressible(mfbac, mfbbc, mfbcc, cy, cy2, c9o2 , c2o9 , oneMinusRho); - backwardInverseChimeraWithKincompressible(mfcaa, mfcba, mfcca, cy, cy2, c6o1 , c1o6 , oneMinusRho); - backwardChimera( mfcab, mfcbb, mfccb, cy, cy2); - backwardInverseChimeraWithKincompressible(mfcac, mfcbc, mfccc, cy, cy2, c18o1, c1o18, oneMinusRho); - - //////////////////////////////////////////////////////////////////////////////////// - // Z - Dir - backwardInverseChimeraWithKincompressible(mfaaa, mfaab, mfaac, cz, cz2, c36o1, c1o36, oneMinusRho); - backwardInverseChimeraWithKincompressible(mfaba, mfabb, mfabc, cz, cz2, c9o1 , c1o9 , oneMinusRho); - backwardInverseChimeraWithKincompressible(mfaca, mfacb, mfacc, cz, cz2, c36o1, c1o36, oneMinusRho); - backwardInverseChimeraWithKincompressible(mfbaa, mfbab, mfbac, cz, cz2, c9o1 , c1o9 , oneMinusRho); - backwardInverseChimeraWithKincompressible(mfbba, mfbbb, mfbbc, cz, cz2, c9o4 , c4o9 , oneMinusRho); - backwardInverseChimeraWithKincompressible(mfbca, mfbcb, mfbcc, cz, cz2, c9o1 , c1o9 , oneMinusRho); - backwardInverseChimeraWithKincompressible(mfcaa, mfcab, mfcac, cz, cz2, c36o1, c1o36, oneMinusRho); - backwardInverseChimeraWithKincompressible(mfcba, mfcbb, mfcbc, cz, cz2, c9o1 , c1o9 , oneMinusRho); - backwardInverseChimeraWithKincompressible(mfcca, mfccb, mfccc, cz, cz2, c36o1, c1o36, oneMinusRho); - - //////////////////////////////////////////////////////////////////////////////////// - //! - Write distributions: style of reading and writing the distributions from/to - //! stored arrays dependent on timestep is based on the esoteric twist algorithm - //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a> - //! - (distAD.f[dirE])[k] = mfabb; - (distAD.f[dirW])[kw] = mfcbb; - (distAD.f[dirN])[k] = mfbab; - (distAD.f[dirS])[ks] = mfbcb; - (distAD.f[dirT])[k] = mfbba; - (distAD.f[dirB])[kb] = mfbbc; - (distAD.f[dirNE])[k] = mfaab; - (distAD.f[dirSW])[ksw] = mfccb; - (distAD.f[dirSE])[ks] = mfacb; - (distAD.f[dirNW])[kw] = mfcab; - (distAD.f[dirTE])[k] = mfaba; - (distAD.f[dirBW])[kbw] = mfcbc; - (distAD.f[dirBE])[kb] = mfabc; - (distAD.f[dirTW])[kw] = mfcba; - (distAD.f[dirTN])[k] = mfbaa; - (distAD.f[dirBS])[kbs] = mfbcc; - (distAD.f[dirBN])[kb] = mfbac; - (distAD.f[dirTS])[ks] = mfbca; - (distAD.f[dirREST])[k] = mfbbb; - (distAD.f[dirTNE])[k] = mfaaa; - (distAD.f[dirTSE])[ks] = mfaca; - (distAD.f[dirBNE])[kb] = mfaac; - (distAD.f[dirBSE])[kbs] = mfacc; - (distAD.f[dirTNW])[kw] = mfcaa; - (distAD.f[dirTSW])[ksw] = mfcca; - (distAD.f[dirBNW])[kbw] = mfcac; - (distAD.f[dirBSW])[kbsw] = mfccc; - } -} -//////////////////////////////////////////////////////////////////////////////// diff --git a/src/gpu/VirtualFluids_GPU/GPU/GPU_Interface.h b/src/gpu/VirtualFluids_GPU/GPU/GPU_Interface.h index d8f28aa6b1a4502894720151f826226abe60d936..bbe25173a87369a1dbf082f77005aa8076e7ebd3 100644 --- a/src/gpu/VirtualFluids_GPU/GPU/GPU_Interface.h +++ b/src/gpu/VirtualFluids_GPU/GPU/GPU_Interface.h @@ -133,21 +133,6 @@ extern "C" void FactorizedCentralMomentsAdvectionDiffusionDeviceKernel( real* forces, bool isEvenTimestep); -////////////////////////////////////////////////////////////////////////// -//! \brief Advection Diffusion kernel -extern "C" void CentralMomentsAdvectionDiffusionDeviceKernel( - uint numberOfThreads, - real omegaD, - uint* typeOfGridNode, - uint* neighborX, - uint* neighborY, - uint* neighborZ, - real* distributions, - real* distributionsAD, - int size_Mat, - real* forces, - bool isEvenTimestep); - ////////////////////////////////////////////////////////////////////////// //! \brief initialize the Advection Diffusion distribution functions extern "C" void InitADdevice( diff --git a/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh b/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh index 92ad82499513e61338a39b40150806c3a383ce1e..fde257f57fb30615829042f2eccfe0ca6dcc62cf 100644 --- a/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh +++ b/src/gpu/VirtualFluids_GPU/GPU/GPU_Kernels.cuh @@ -129,20 +129,6 @@ extern "C" __global__ void Factorized_Central_Moments_Advection_Diffusion_Device real* forces, bool isEvenTimestep); -////////////////////////////////////////////////////////////////////////// -//! \brief \ref Advection_Diffusion_Device_Kernel : central moments for Advection Diffusion Equation -extern "C" __global__ void Central_Moments_Advection_Diffusion_Device_Kernel( - real omegaD, - uint* typeOfGridNode, - uint* neighborX, - uint* neighborY, - uint* neighborZ, - real* distributions, - real* distributionsAD, - int size_Mat, - real* forces, - bool isEvenTimestep); - ////////////////////////////////////////////////////////////////////////// //! \brief \ref InitAD : device function to initialize the distributionsAD extern "C" __global__ void InitAD( diff --git a/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu b/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu index 834ac1f4aa5697bcb35d25f82d549168b70381ce..3791556fb7e9d957f0771ab746a0d002e7a92c94 100644 --- a/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu +++ b/src/gpu/VirtualFluids_GPU/GPU/LBMKernel.cu @@ -244,37 +244,6 @@ extern "C" void FactorizedCentralMomentsAdvectionDiffusionDeviceKernel( getLastCudaError("Factorized_Central_Moments_Advection_Diffusion_Device_Kernel execution failed"); } ////////////////////////////////////////////////////////////////////////// -extern "C" void CentralMomentsAdvectionDiffusionDeviceKernel( - uint numberOfThreads, - real omegaD, - uint* typeOfGridNode, - uint* neighborX, - uint* neighborY, - uint* neighborZ, - real* distributions, - real* distributionsAD, - int size_Mat, - real* forces, - bool isEvenTimestep) -{ - int Grid = (size_Mat / numberOfThreads) + 1; - dim3 grid(Grid, 1, 1); - dim3 threads(numberOfThreads, 1, 1); - - Central_Moments_Advection_Diffusion_Device_Kernel <<< grid, threads >>> ( - omegaD, - typeOfGridNode, - neighborX, - neighborY, - neighborZ, - distributions, - distributionsAD, - size_Mat, - forces, - isEvenTimestep); - getLastCudaError("Central_Moments_Advection_Diffusion_Device_Kernel execution failed"); -} -////////////////////////////////////////////////////////////////////////// extern "C" void InitADdevice( uint numberOfThreads, uint* neighborX, diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp index a32bbd70e396dae3dc5699690e4b5d90c621a6dc..14ed3fdbb63e0061b835d5f9235a442366dccc8f 100644 --- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp +++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp @@ -341,12 +341,12 @@ std::function<void(real, real, real, real&, real&, real&, real&)>& Parameter::ge } // initial condition concentration -void Parameter::setInitialConditionAD(std::function<void(real, real, real, real&, real)> initialConditionAD) +void Parameter::setInitialConditionAD(std::function<void(real, real, real, real&)> initialConditionAD) { this->initialConditionAD = initialConditionAD; } -std::function<void(real, real, real, real&, real)>& Parameter::getInitialConditionAD() +std::function<void(real, real, real, real&)>& Parameter::getInitialConditionAD() { return this->initialConditionAD; }