diff options
Diffstat (limited to 'extern')
-rw-r--r-- | extern/mantaflow/CMakeLists.txt | 1 | ||||
-rw-r--r-- | extern/mantaflow/helper/util/rcmatrix.h | 2 | ||||
-rw-r--r-- | extern/mantaflow/preprocessed/conjugategrad.cpp | 94 | ||||
-rw-r--r-- | extern/mantaflow/preprocessed/conjugategrad.h | 466 | ||||
-rw-r--r-- | extern/mantaflow/preprocessed/general.h | 35 | ||||
-rw-r--r-- | extern/mantaflow/preprocessed/gitinfo.h | 2 | ||||
-rw-r--r-- | extern/mantaflow/preprocessed/plugin/extforces.cpp | 25 | ||||
-rw-r--r-- | extern/mantaflow/preprocessed/plugin/pressure.cpp | 14 | ||||
-rw-r--r-- | extern/mantaflow/preprocessed/plugin/viscosity.cpp | 1430 | ||||
-rw-r--r-- | extern/mantaflow/preprocessed/plugin/vortexplugins.cpp | 4 | ||||
-rw-r--r-- | extern/mantaflow/preprocessed/plugin/waves.cpp | 13 | ||||
-rw-r--r-- | extern/mantaflow/preprocessed/registration.cpp | 2 |
12 files changed, 1965 insertions, 123 deletions
diff --git a/extern/mantaflow/CMakeLists.txt b/extern/mantaflow/CMakeLists.txt index ccf272650e3..82bf95a9742 100644 --- a/extern/mantaflow/CMakeLists.txt +++ b/extern/mantaflow/CMakeLists.txt @@ -200,6 +200,7 @@ set(SRC ${MANTA_PP}/plugin/ptsplugins.cpp ${MANTA_PP}/plugin/secondaryparticles.cpp ${MANTA_PP}/plugin/surfaceturbulence.cpp + ${MANTA_PP}/plugin/viscosity.cpp ${MANTA_PP}/plugin/vortexplugins.cpp ${MANTA_PP}/plugin/waveletturbulence.cpp ${MANTA_PP}/plugin/waves.cpp diff --git a/extern/mantaflow/helper/util/rcmatrix.h b/extern/mantaflow/helper/util/rcmatrix.h index f1f0efe6416..330fd1f64f7 100644 --- a/extern/mantaflow/helper/util/rcmatrix.h +++ b/extern/mantaflow/helper/util/rcmatrix.h @@ -1035,7 +1035,7 @@ template<class N, class T> struct RCFixedMatrix { typedef RCMatrix<int, Real> Matrix; typedef RCFixedMatrix<int, Real> FixedMatrix; -} // namespace Manta +} #undef parallel_for #undef parallel_end diff --git a/extern/mantaflow/preprocessed/conjugategrad.cpp b/extern/mantaflow/preprocessed/conjugategrad.cpp index da82a412a97..bdcceb29520 100644 --- a/extern/mantaflow/preprocessed/conjugategrad.cpp +++ b/extern/mantaflow/preprocessed/conjugategrad.cpp @@ -397,7 +397,7 @@ struct UpdateSearchVec : public KernelBase { }; //***************************************************************************** -// CG class +// CG class template<class APPLYMAT> GridCg<APPLYMAT>::GridCg(Grid<Real> &dst, @@ -406,10 +406,8 @@ GridCg<APPLYMAT>::GridCg(Grid<Real> &dst, Grid<Real> &search, const FlagGrid &flags, Grid<Real> &tmp, - Grid<Real> *pA0, - Grid<Real> *pAi, - Grid<Real> *pAj, - Grid<Real> *pAk) + std::vector<Grid<Real> *> matrixAVec, + std::vector<Grid<Real> *> rhsVec) : GridCgInterface(), mInited(false), mIterations(0), @@ -419,10 +417,8 @@ GridCg<APPLYMAT>::GridCg(Grid<Real> &dst, mSearch(search), mFlags(flags), mTmp(tmp), - mpA0(pA0), - mpAi(pAi), - mpAj(pAj), - mpAk(pAk), + mMatrixA(matrixAVec), + mVecRhs(rhsVec), mPcMethod(PC_None), mpPCA0(nullptr), mpPCAi(nullptr), @@ -445,19 +441,37 @@ template<class APPLYMAT> void GridCg<APPLYMAT>::doInit() if (mPcMethod == PC_ICP) { assertMsg(mDst.is3D(), "ICP only supports 3D grids so far"); - InitPreconditionIncompCholesky( - mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk); - ApplyPreconditionIncompCholesky( - mTmp, mResidual, mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk); + InitPreconditionIncompCholesky(mFlags, + *mpPCA0, + *mpPCAi, + *mpPCAj, + *mpPCAk, + *mMatrixA[0], + *mMatrixA[1], + *mMatrixA[2], + *mMatrixA[3]); + ApplyPreconditionIncompCholesky(mTmp, + mResidual, + mFlags, + *mpPCA0, + *mpPCAi, + *mpPCAj, + *mpPCAk, + *mMatrixA[0], + *mMatrixA[1], + *mMatrixA[2], + *mMatrixA[3]); } else if (mPcMethod == PC_mICP) { assertMsg(mDst.is3D(), "mICP only supports 3D grids so far"); - InitPreconditionModifiedIncompCholesky2(mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk); + InitPreconditionModifiedIncompCholesky2( + mFlags, *mpPCA0, *mMatrixA[0], *mMatrixA[1], *mMatrixA[2], *mMatrixA[3]); ApplyPreconditionModifiedIncompCholesky2( - mTmp, mResidual, mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk); + mTmp, mResidual, mFlags, *mpPCA0, *mMatrixA[0], *mMatrixA[1], *mMatrixA[2], *mMatrixA[3]); } else if (mPcMethod == PC_MGP) { - InitPreconditionMultigrid(mMG, *mpA0, *mpAi, *mpAj, *mpAk, mAccuracy); + InitPreconditionMultigrid( + mMG, *mMatrixA[0], *mMatrixA[1], *mMatrixA[2], *mMatrixA[3], mAccuracy); ApplyPreconditionMultigrid(mMG, mTmp, mResidual); } else { @@ -465,7 +479,6 @@ template<class APPLYMAT> void GridCg<APPLYMAT>::doInit() } mSearch.copyFrom(mTmp); - mSigma = GridDotProduct(mTmp, mResidual); } @@ -480,7 +493,7 @@ template<class APPLYMAT> bool GridCg<APPLYMAT>::iterate() // this could reinterpret the mpA pointers (not so clean right now) // tmp = applyMat(search) - APPLYMAT(mFlags, mTmp, mSearch, *mpA0, *mpAi, *mpAj, *mpAk); + APPLYMAT(mFlags, mTmp, mSearch, mMatrixA, mVecRhs); // alpha = sigma/dot(tmp, search) Real dp = GridDotProduct(mTmp, mSearch); @@ -492,11 +505,20 @@ template<class APPLYMAT> bool GridCg<APPLYMAT>::iterate() gridScaledAdd<Real, Real>(mResidual, mTmp, -alpha); // residual += tmp * -alpha if (mPcMethod == PC_ICP) - ApplyPreconditionIncompCholesky( - mTmp, mResidual, mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk); + ApplyPreconditionIncompCholesky(mTmp, + mResidual, + mFlags, + *mpPCA0, + *mpPCAi, + *mpPCAj, + *mpPCAk, + *mMatrixA[0], + *mMatrixA[1], + *mMatrixA[2], + *mMatrixA[3]); else if (mPcMethod == PC_mICP) ApplyPreconditionModifiedIncompCholesky2( - mTmp, mResidual, mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk); + mTmp, mResidual, mFlags, *mpPCA0, *mMatrixA[0], *mMatrixA[1], *mMatrixA[2], *mMatrixA[3]); else if (mPcMethod == PC_MGP) ApplyPreconditionMultigrid(mMG, mTmp, mResidual); else @@ -584,13 +606,15 @@ void GridCg<APPLYMAT>::setMGPreconditioner(PreconditionType method, GridMg *MG) assertMsg(method == PC_MGP, "GridCg<APPLYMAT>::setMGPreconditioner: Invalid method specified."); mPcMethod = method; - mMG = MG; } // explicit instantiation template class GridCg<ApplyMatrix>; template class GridCg<ApplyMatrix2D>; +template class GridCg<ApplyMatrixViscosityU>; +template class GridCg<ApplyMatrixViscosityV>; +template class GridCg<ApplyMatrixViscosityW>; //***************************************************************************** // diffusion for real and vec grids, e.g. for viscosity @@ -638,10 +662,15 @@ void cgSolveDiffusion(const FlagGrid &flags, if (grid.getType() & GridBase::TypeReal) { Grid<Real> &u = ((Grid<Real> &)grid); rhs.copyFrom(u); - if (flags.is3D()) - gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak); - else - gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak); + vector<Grid<Real> *> matA{&A0, &Ai, &Aj}; + + if (flags.is3D()) { + matA.push_back(&Ak); + gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, matA); + } + else { + gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, matA); + } gcg->setAccuracy(cgAccuracy); gcg->solve(maxIter); @@ -653,12 +682,17 @@ void cgSolveDiffusion(const FlagGrid &flags, else if ((grid.getType() & GridBase::TypeVec3) || (grid.getType() & GridBase::TypeMAC)) { Grid<Vec3> &vec = ((Grid<Vec3> &)grid); Grid<Real> u(parent); + vector<Grid<Real> *> matA{&A0, &Ai, &Aj}; // core solve is same as for a regular real grid - if (flags.is3D()) - gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak); - else - gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak); + if (flags.is3D()) { + matA.push_back(&Ak); + gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, matA); + } + else { + gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, matA); + } + gcg->setAccuracy(cgAccuracy); // diffuse every component separately diff --git a/extern/mantaflow/preprocessed/conjugategrad.h b/extern/mantaflow/preprocessed/conjugategrad.h index 58ccff28179..4974aa6a4d6 100644 --- a/extern/mantaflow/preprocessed/conjugategrad.h +++ b/extern/mantaflow/preprocessed/conjugategrad.h @@ -78,13 +78,9 @@ template<class APPLYMAT> class GridCg : public GridCgInterface { Grid<Real> &search, const FlagGrid &flags, Grid<Real> &tmp, - Grid<Real> *A0, - Grid<Real> *pAi, - Grid<Real> *pAj, - Grid<Real> *pAk); - ~GridCg() - { - } + std::vector<Grid<Real> *> matrixAVec, + std::vector<Grid<Real> *> rhsVec = {}); + ~GridCg(){}; void doInit(); bool iterate(); @@ -133,7 +129,10 @@ template<class APPLYMAT> class GridCg : public GridCgInterface { const FlagGrid &mFlags; Grid<Real> &mTmp; - Grid<Real> *mpA0, *mpAi, *mpAj, *mpAk; + //! shape of A matrix defined here (e.g. diagonal, positive neighbor cells, etc) + std::vector<Grid<Real> *> mMatrixA; + //! shape of rhs vector defined here (e.g. 1 rhs for regular fluids solve, 3 rhs for viscosity) + std::vector<Grid<Real> *> mVecRhs; PreconditionType mPcMethod; //! preconditioning grids @@ -154,11 +153,9 @@ struct ApplyMatrix : public KernelBase { ApplyMatrix(const FlagGrid &flags, Grid<Real> &dst, const Grid<Real> &src, - Grid<Real> &A0, - Grid<Real> &Ai, - Grid<Real> &Aj, - Grid<Real> &Ak) - : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), A0(A0), Ai(Ai), Aj(Aj), Ak(Ak) + const std::vector<Grid<Real> *> matrixA, + const std::vector<Grid<Real> *> vecRhs) + : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs) { runMessage(); run(); @@ -167,11 +164,18 @@ struct ApplyMatrix : public KernelBase { const FlagGrid &flags, Grid<Real> &dst, const Grid<Real> &src, - Grid<Real> &A0, - Grid<Real> &Ai, - Grid<Real> &Aj, - Grid<Real> &Ak) const + const std::vector<Grid<Real> *> matrixA, + const std::vector<Grid<Real> *> vecRhs) const { + unusedParameter(vecRhs); // Not needed in this matrix application + + if (matrixA.size() != 4) + errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step"); + Grid<Real> &A0 = *matrixA[0]; + Grid<Real> &Ai = *matrixA[1]; + Grid<Real> &Aj = *matrixA[2]; + Grid<Real> &Ak = *matrixA[3]; + if (!flags.isFluid(idx)) { dst[idx] = src[idx]; return; @@ -196,26 +200,16 @@ struct ApplyMatrix : public KernelBase { return src; } typedef Grid<Real> type2; - inline Grid<Real> &getArg3() - { - return A0; - } - typedef Grid<Real> type3; - inline Grid<Real> &getArg4() - { - return Ai; - } - typedef Grid<Real> type4; - inline Grid<Real> &getArg5() + inline const std::vector<Grid<Real> *> &getArg3() { - return Aj; + return matrixA; } - typedef Grid<Real> type5; - inline Grid<Real> &getArg6() + typedef std::vector<Grid<Real> *> type3; + inline const std::vector<Grid<Real> *> &getArg4() { - return Ak; + return vecRhs; } - typedef Grid<Real> type6; + typedef std::vector<Grid<Real> *> type4; void runMessage() { debMsg("Executing kernel ApplyMatrix ", 3); @@ -226,7 +220,7 @@ struct ApplyMatrix : public KernelBase { void operator()(const tbb::blocked_range<IndexInt> &__r) const { for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++) - op(idx, flags, dst, src, A0, Ai, Aj, Ak); + op(idx, flags, dst, src, matrixA, vecRhs); } void run() { @@ -235,10 +229,8 @@ struct ApplyMatrix : public KernelBase { const FlagGrid &flags; Grid<Real> &dst; const Grid<Real> &src; - Grid<Real> &A0; - Grid<Real> &Ai; - Grid<Real> &Aj; - Grid<Real> &Ak; + const std::vector<Grid<Real> *> matrixA; + const std::vector<Grid<Real> *> vecRhs; }; //! Kernel: Apply symmetric stored Matrix. 2D version @@ -247,11 +239,9 @@ struct ApplyMatrix2D : public KernelBase { ApplyMatrix2D(const FlagGrid &flags, Grid<Real> &dst, const Grid<Real> &src, - Grid<Real> &A0, - Grid<Real> &Ai, - Grid<Real> &Aj, - Grid<Real> &Ak) - : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), A0(A0), Ai(Ai), Aj(Aj), Ak(Ak) + const std::vector<Grid<Real> *> matrixA, + const std::vector<Grid<Real> *> vecRhs) + : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs) { runMessage(); run(); @@ -260,12 +250,16 @@ struct ApplyMatrix2D : public KernelBase { const FlagGrid &flags, Grid<Real> &dst, const Grid<Real> &src, - Grid<Real> &A0, - Grid<Real> &Ai, - Grid<Real> &Aj, - Grid<Real> &Ak) const + const std::vector<Grid<Real> *> matrixA, + const std::vector<Grid<Real> *> vecRhs) const { - unusedParameter(Ak); // only there for parameter compatibility with ApplyMatrix + unusedParameter(vecRhs); // Not needed in this matrix application + + if (matrixA.size() != 3) + errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step"); + Grid<Real> &A0 = *matrixA[0]; + Grid<Real> &Ai = *matrixA[1]; + Grid<Real> &Aj = *matrixA[2]; if (!flags.isFluid(idx)) { dst[idx] = src[idx]; @@ -290,51 +284,387 @@ struct ApplyMatrix2D : public KernelBase { return src; } typedef Grid<Real> type2; - inline Grid<Real> &getArg3() + inline const std::vector<Grid<Real> *> &getArg3() { - return A0; + return matrixA; } - typedef Grid<Real> type3; - inline Grid<Real> &getArg4() + typedef std::vector<Grid<Real> *> type3; + inline const std::vector<Grid<Real> *> &getArg4() { - return Ai; + return vecRhs; } - typedef Grid<Real> type4; - inline Grid<Real> &getArg5() + typedef std::vector<Grid<Real> *> type4; + void runMessage() { - return Aj; + debMsg("Executing kernel ApplyMatrix2D ", 3); + debMsg("Kernel range" + << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ", + 4); + }; + void operator()(const tbb::blocked_range<IndexInt> &__r) const + { + for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++) + op(idx, flags, dst, src, matrixA, vecRhs); } - typedef Grid<Real> type5; - inline Grid<Real> &getArg6() + void run() { - return Ak; + tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this); + } + const FlagGrid &flags; + Grid<Real> &dst; + const Grid<Real> &src; + const std::vector<Grid<Real> *> matrixA; + const std::vector<Grid<Real> *> vecRhs; +}; + +struct ApplyMatrixViscosityU : public KernelBase { + ApplyMatrixViscosityU(const FlagGrid &flags, + Grid<Real> &dst, + const Grid<Real> &src, + const std::vector<Grid<Real> *> matrixA, + const std::vector<Grid<Real> *> vecRhs) + : KernelBase(&flags, 1), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs) + { + runMessage(); + run(); } - typedef Grid<Real> type6; + inline void op(int i, + int j, + int k, + const FlagGrid &flags, + Grid<Real> &dst, + const Grid<Real> &src, + const std::vector<Grid<Real> *> matrixA, + const std::vector<Grid<Real> *> vecRhs) const + { + if (matrixA.size() != 15) + errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step"); + Grid<Real> &A0 = *matrixA[0]; + Grid<Real> &Aplusi = *matrixA[1]; + Grid<Real> &Aplusj = *matrixA[2]; + Grid<Real> &Aplusk = *matrixA[3]; + Grid<Real> &Aminusi = *matrixA[4]; + Grid<Real> &Aminusj = *matrixA[5]; + Grid<Real> &Aminusk = *matrixA[6]; + + if (vecRhs.size() != 2) + errMsg("ConjugatedGrad: Invalid rhs vector in apply matrix step"); + Grid<Real> &srcV = *vecRhs[0]; + Grid<Real> &srcW = *vecRhs[1]; + + dst(i, j, k) = src(i, j, k) * A0(i, j, k) + src(i + 1, j, k) * Aplusi(i, j, k) + + src(i, j + 1, k) * Aplusj(i, j, k) + src(i, j, k + 1) * Aplusk(i, j, k) + + src(i - 1, j, k) * Aminusi(i, j, k) + src(i, j - 1, k) * Aminusj(i, j, k) + + src(i, j, k - 1) * Aminusk(i, j, k); + + dst(i, j, k) += srcV(i, j + 1, k) * (*matrixA[7])(i, j, k) + + srcV(i - 1, j + 1, k) * (*matrixA[8])(i, j, k) + + srcV(i, j, k) * (*matrixA[9])(i, j, k) + + srcV(i - 1, j, k) * (*matrixA[10])(i, j, k) + + srcW(i, j, k + 1) * (*matrixA[11])(i, j, k) + + srcW(i - 1, j, k + 1) * (*matrixA[12])(i, j, k) + + srcW(i, j, k) * (*matrixA[13])(i, j, k) + + srcW(i - 1, j, k) * (*matrixA[14])(i, j, k); + } + inline const FlagGrid &getArg0() + { + return flags; + } + typedef FlagGrid type0; + inline Grid<Real> &getArg1() + { + return dst; + } + typedef Grid<Real> type1; + inline const Grid<Real> &getArg2() + { + return src; + } + typedef Grid<Real> type2; + inline const std::vector<Grid<Real> *> &getArg3() + { + return matrixA; + } + typedef std::vector<Grid<Real> *> type3; + inline const std::vector<Grid<Real> *> &getArg4() + { + return vecRhs; + } + typedef std::vector<Grid<Real> *> type4; void runMessage() { - debMsg("Executing kernel ApplyMatrix2D ", 3); + debMsg("Executing kernel ApplyMatrixViscosityU ", 3); debMsg("Kernel range" << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ", 4); }; void operator()(const tbb::blocked_range<IndexInt> &__r) const { - for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++) - op(idx, flags, dst, src, A0, Ai, Aj, Ak); + const int _maxX = maxX; + const int _maxY = maxY; + if (maxZ > 1) { + for (int k = __r.begin(); k != (int)__r.end(); k++) + for (int j = 1; j < _maxY; j++) + for (int i = 1; i < _maxX; i++) + op(i, j, k, flags, dst, src, matrixA, vecRhs); + } + else { + const int k = 0; + for (int j = __r.begin(); j != (int)__r.end(); j++) + for (int i = 1; i < _maxX; i++) + op(i, j, k, flags, dst, src, matrixA, vecRhs); + } } void run() { - tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this); + if (maxZ > 1) + tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this); + else + tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this); } const FlagGrid &flags; Grid<Real> &dst; const Grid<Real> &src; - Grid<Real> &A0; - Grid<Real> &Ai; - Grid<Real> &Aj; - Grid<Real> &Ak; + const std::vector<Grid<Real> *> matrixA; + const std::vector<Grid<Real> *> vecRhs; }; +struct ApplyMatrixViscosityV : public KernelBase { + ApplyMatrixViscosityV(const FlagGrid &flags, + Grid<Real> &dst, + const Grid<Real> &src, + const std::vector<Grid<Real> *> matrixA, + const std::vector<Grid<Real> *> vecRhs) + : KernelBase(&flags, 1), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs) + { + runMessage(); + run(); + } + inline void op(int i, + int j, + int k, + const FlagGrid &flags, + Grid<Real> &dst, + const Grid<Real> &src, + const std::vector<Grid<Real> *> matrixA, + const std::vector<Grid<Real> *> vecRhs) const + { + if (matrixA.size() != 15) + errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step"); + Grid<Real> &A0 = *matrixA[0]; + Grid<Real> &Aplusi = *matrixA[1]; + Grid<Real> &Aplusj = *matrixA[2]; + Grid<Real> &Aplusk = *matrixA[3]; + Grid<Real> &Aminusi = *matrixA[4]; + Grid<Real> &Aminusj = *matrixA[5]; + Grid<Real> &Aminusk = *matrixA[6]; + + if (vecRhs.size() != 2) + errMsg("ConjugatedGrad: Invalid rhs vector in apply matrix step"); + Grid<Real> &srcU = *vecRhs[0]; + Grid<Real> &srcW = *vecRhs[1]; + + dst(i, j, k) = src(i, j, k) * A0(i, j, k) + src(i + 1, j, k) * Aplusi(i, j, k) + + src(i, j + 1, k) * Aplusj(i, j, k) + src(i, j, k + 1) * Aplusk(i, j, k) + + src(i - 1, j, k) * Aminusi(i, j, k) + src(i, j - 1, k) * Aminusj(i, j, k) + + src(i, j, k - 1) * Aminusk(i, j, k); + + dst(i, j, k) += srcU(i + 1, j, k) * (*matrixA[7])(i, j, k) + + srcU(i + 1, j - 1, k) * (*matrixA[8])(i, j, k) + + srcU(i, j, k) * (*matrixA[9])(i, j, k) + + srcU(i, j - 1, k) * (*matrixA[10])(i, j, k) + + srcW(i, j, k + 1) * (*matrixA[11])(i, j, k) + + srcW(i, j - 1, k + 1) * (*matrixA[12])(i, j, k) + + srcW(i, j, k) * (*matrixA[13])(i, j, k) + + srcW(i, j - 1, k) * (*matrixA[14])(i, j, k); + } + inline const FlagGrid &getArg0() + { + return flags; + } + typedef FlagGrid type0; + inline Grid<Real> &getArg1() + { + return dst; + } + typedef Grid<Real> type1; + inline const Grid<Real> &getArg2() + { + return src; + } + typedef Grid<Real> type2; + inline const std::vector<Grid<Real> *> &getArg3() + { + return matrixA; + } + typedef std::vector<Grid<Real> *> type3; + inline const std::vector<Grid<Real> *> &getArg4() + { + return vecRhs; + } + typedef std::vector<Grid<Real> *> type4; + void runMessage() + { + debMsg("Executing kernel ApplyMatrixViscosityV ", 3); + debMsg("Kernel range" + << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ", + 4); + }; + void operator()(const tbb::blocked_range<IndexInt> &__r) const + { + const int _maxX = maxX; + const int _maxY = maxY; + if (maxZ > 1) { + for (int k = __r.begin(); k != (int)__r.end(); k++) + for (int j = 1; j < _maxY; j++) + for (int i = 1; i < _maxX; i++) + op(i, j, k, flags, dst, src, matrixA, vecRhs); + } + else { + const int k = 0; + for (int j = __r.begin(); j != (int)__r.end(); j++) + for (int i = 1; i < _maxX; i++) + op(i, j, k, flags, dst, src, matrixA, vecRhs); + } + } + void run() + { + if (maxZ > 1) + tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this); + else + tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this); + } + const FlagGrid &flags; + Grid<Real> &dst; + const Grid<Real> &src; + const std::vector<Grid<Real> *> matrixA; + const std::vector<Grid<Real> *> vecRhs; +}; + +struct ApplyMatrixViscosityW : public KernelBase { + ApplyMatrixViscosityW(const FlagGrid &flags, + Grid<Real> &dst, + const Grid<Real> &src, + const std::vector<Grid<Real> *> matrixA, + const std::vector<Grid<Real> *> vecRhs) + : KernelBase(&flags, 1), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs) + { + runMessage(); + run(); + } + inline void op(int i, + int j, + int k, + const FlagGrid &flags, + Grid<Real> &dst, + const Grid<Real> &src, + const std::vector<Grid<Real> *> matrixA, + const std::vector<Grid<Real> *> vecRhs) const + { + if (matrixA.size() != 15) + errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step"); + Grid<Real> &A0 = *matrixA[0]; + Grid<Real> &Aplusi = *matrixA[1]; + Grid<Real> &Aplusj = *matrixA[2]; + Grid<Real> &Aplusk = *matrixA[3]; + Grid<Real> &Aminusi = *matrixA[4]; + Grid<Real> &Aminusj = *matrixA[5]; + Grid<Real> &Aminusk = *matrixA[6]; + + if (vecRhs.size() != 2) + errMsg("ConjugatedGrad: Invalid rhs vector in apply matrix step"); + Grid<Real> &srcU = *vecRhs[0]; + Grid<Real> &srcV = *vecRhs[1]; + + dst(i, j, k) = src(i, j, k) * A0(i, j, k) + src(i + 1, j, k) * Aplusi(i, j, k) + + src(i, j + 1, k) * Aplusj(i, j, k) + src(i, j, k + 1) * Aplusk(i, j, k) + + src(i - 1, j, k) * Aminusi(i, j, k) + src(i, j - 1, k) * Aminusj(i, j, k) + + src(i, j, k - 1) * Aminusk(i, j, k); + + dst(i, j, k) += srcU(i + 1, j, k) * (*matrixA[7])(i, j, k) + + srcU(i + 1, j, k - 1) * (*matrixA[8])(i, j, k) + + srcU(i, j, k) * (*matrixA[9])(i, j, k) + + srcU(i, j, k - 1) * (*matrixA[10])(i, j, k) + + srcV(i, j + 1, k) * (*matrixA[11])(i, j, k) + + srcV(i, j + 1, k - 1) * (*matrixA[12])(i, j, k) + + srcV(i, j, k) * (*matrixA[13])(i, j, k) + + srcV(i, j, k - 1) * (*matrixA[14])(i, j, k); + } + inline const FlagGrid &getArg0() + { + return flags; + } + typedef FlagGrid type0; + inline Grid<Real> &getArg1() + { + return dst; + } + typedef Grid<Real> type1; + inline const Grid<Real> &getArg2() + { + return src; + } + typedef Grid<Real> type2; + inline const std::vector<Grid<Real> *> &getArg3() + { + return matrixA; + } + typedef std::vector<Grid<Real> *> type3; + inline const std::vector<Grid<Real> *> &getArg4() + { + return vecRhs; + } + typedef std::vector<Grid<Real> *> type4; + void runMessage() + { + debMsg("Executing kernel ApplyMatrixViscosityW ", 3); + debMsg("Kernel range" + << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ", + 4); + }; + void operator()(const tbb::blocked_range<IndexInt> &__r) const + { + const int _maxX = maxX; + const int _maxY = maxY; + if (maxZ > 1) { + for (int k = __r.begin(); k != (int)__r.end(); k++) + for (int j = 1; j < _maxY; j++) + for (int i = 1; i < _maxX; i++) + op(i, j, k, flags, dst, src, matrixA, vecRhs); + } + else { + const int k = 0; + for (int j = __r.begin(); j != (int)__r.end(); j++) + for (int i = 1; i < _maxX; i++) + op(i, j, k, flags, dst, src, matrixA, vecRhs); + } + } + void run() + { + if (maxZ > 1) + tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this); + else + tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this); + } + const FlagGrid &flags; + Grid<Real> &dst; + const Grid<Real> &src; + const std::vector<Grid<Real> *> matrixA; + const std::vector<Grid<Real> *> vecRhs; +}; + +/* NOTE: Use this template for new matrix application kernels + +//! Template for matrix application kernels +KERNEL() +void ApplyMatrixTemplate (const FlagGrid& flags, Grid<Real>& dst, const Grid<Real>& src, + const std::vector<Grid<Real> *> matrixA, const std::vector<Grid<Real> *> vecRhs) +{ + // The kernel must define how to use the grids from the matrixA and vecRhs lists +} + +*/ + //! Kernel: Construct the matrix for the poisson equation struct MakeLaplaceMatrix : public KernelBase { diff --git a/extern/mantaflow/preprocessed/general.h b/extern/mantaflow/preprocessed/general.h index 7a840517cef..50eac71e87e 100644 --- a/extern/mantaflow/preprocessed/general.h +++ b/extern/mantaflow/preprocessed/general.h @@ -42,7 +42,7 @@ inline void updateQtGui(bool full, int frame, float time, const std::string &cur # ifdef _DEBUG # define DEBUG 1 # endif // _DEBUG -#endif // DEBUG +#endif // DEBUG // Standard exception class Error : public std::exception { @@ -242,6 +242,39 @@ inline bool c_isnan(float c) return d != d; } +//! Swap so that a<b +template<class T> inline void sort(T &a, T &b) +{ + if (a > b) + std::swap(a, b); +} + +//! Swap so that a<b<c +template<class T> inline void sort(T &a, T &b, T &c) +{ + if (a > b) + std::swap(a, b); + if (a > c) + std::swap(a, c); + if (b > c) + std::swap(b, c); +} + +//! Swap so that a<b<c<d +template<class T> inline void sort(T &a, T &b, T &c, T &d) +{ + if (a > b) + std::swap(a, b); + if (c > d) + std::swap(c, d); + if (a > c) + std::swap(a, c); + if (b > d) + std::swap(b, d); + if (b > c) + std::swap(b, c); +} + } // namespace Manta #endif diff --git a/extern/mantaflow/preprocessed/gitinfo.h b/extern/mantaflow/preprocessed/gitinfo.h index a0590e6a0b0..afc34a00d3d 100644 --- a/extern/mantaflow/preprocessed/gitinfo.h +++ b/extern/mantaflow/preprocessed/gitinfo.h @@ -1,3 +1,3 @@ -#define MANTA_GIT_VERSION "commit 327917cd59b03bef3a953b5f58fc1637b3a83e01" +#define MANTA_GIT_VERSION "commit e2285cb9bc492987f728123be6cfc1fe11fe73d6" diff --git a/extern/mantaflow/preprocessed/plugin/extforces.cpp b/extern/mantaflow/preprocessed/plugin/extforces.cpp index 558008afe3b..88935fa7ae9 100644 --- a/extern/mantaflow/preprocessed/plugin/extforces.cpp +++ b/extern/mantaflow/preprocessed/plugin/extforces.cpp @@ -1135,26 +1135,27 @@ struct KnAddForceIfLower : public KernelBase { if (!curFluid && !curEmpty) return; + Real minVal, maxVal, sum; if (flags.isFluid(i - 1, j, k) || (curFluid && flags.isEmpty(i - 1, j, k))) { Real forceMACX = 0.5 * (force(i - 1, j, k).x + force(i, j, k).x); - Real min = std::min(vel(i, j, k).x, forceMACX); - Real max = std::max(vel(i, j, k).x, forceMACX); - Real sum = vel(i, j, k).x + forceMACX; - vel(i, j, k).x = (forceMACX > 0) ? std::min(sum, max) : std::max(sum, min); + minVal = min(vel(i, j, k).x, forceMACX); + maxVal = max(vel(i, j, k).x, forceMACX); + sum = vel(i, j, k).x + forceMACX; + vel(i, j, k).x = (forceMACX > 0) ? min(sum, maxVal) : max(sum, minVal); } if (flags.isFluid(i, j - 1, k) || (curFluid && flags.isEmpty(i, j - 1, k))) { Real forceMACY = 0.5 * (force(i, j - 1, k).y + force(i, j, k).y); - Real min = std::min(vel(i, j, k).y, forceMACY); - Real max = std::max(vel(i, j, k).y, forceMACY); - Real sum = vel(i, j, k).y + forceMACY; - vel(i, j, k).y = (forceMACY > 0) ? std::min(sum, max) : std::max(sum, min); + minVal = min(vel(i, j, k).y, forceMACY); + maxVal = max(vel(i, j, k).y, forceMACY); + sum = vel(i, j, k).y + forceMACY; + vel(i, j, k).y = (forceMACY > 0) ? min(sum, maxVal) : max(sum, minVal); } if (vel.is3D() && (flags.isFluid(i, j, k - 1) || (curFluid && flags.isEmpty(i, j, k - 1)))) { Real forceMACZ = 0.5 * (force(i, j, k - 1).z + force(i, j, k).z); - Real min = std::min(vel(i, j, k).z, forceMACZ); - Real max = std::max(vel(i, j, k).z, forceMACZ); - Real sum = vel(i, j, k).z + forceMACZ; - vel(i, j, k).z = (forceMACZ > 0) ? std::min(sum, max) : std::max(sum, min); + minVal = min(vel(i, j, k).z, forceMACZ); + maxVal = max(vel(i, j, k).z, forceMACZ); + sum = vel(i, j, k).z + forceMACZ; + vel(i, j, k).z = (forceMACZ > 0) ? min(sum, maxVal) : max(sum, minVal); } } inline const FlagGrid &getArg0() diff --git a/extern/mantaflow/preprocessed/plugin/pressure.cpp b/extern/mantaflow/preprocessed/plugin/pressure.cpp index 1100a58db47..593aeb16859 100644 --- a/extern/mantaflow/preprocessed/plugin/pressure.cpp +++ b/extern/mantaflow/preprocessed/plugin/pressure.cpp @@ -1138,11 +1138,15 @@ void solvePressureSystem(Grid<Real> &rhs, // note: the last factor increases the max iterations for 2d, which right now can't use a // preconditioner GridCgInterface *gcg; - if (vel.is3D()) - gcg = new GridCg<ApplyMatrix>(pressure, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak); - else - gcg = new GridCg<ApplyMatrix2D>( - pressure, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak); + vector<Grid<Real> *> matA{&A0, &Ai, &Aj}; + + if (vel.is3D()) { + matA.push_back(&Ak); + gcg = new GridCg<ApplyMatrix>(pressure, rhs, residual, search, flags, tmp, matA); + } + else { + gcg = new GridCg<ApplyMatrix2D>(pressure, rhs, residual, search, flags, tmp, matA); + } gcg->setAccuracy(cgAccuracy); gcg->setUseL2Norm(useL2Norm); diff --git a/extern/mantaflow/preprocessed/plugin/viscosity.cpp b/extern/mantaflow/preprocessed/plugin/viscosity.cpp new file mode 100644 index 00000000000..5dc478b13dd --- /dev/null +++ b/extern/mantaflow/preprocessed/plugin/viscosity.cpp @@ -0,0 +1,1430 @@ + + +// DO NOT EDIT ! +// This file is generated using the MantaFlow preprocessor (prep generate). + +/****************************************************************************** + * + * MantaFlow fluid solver framework + * Copyright 2020 Sebastian Barschkis, Nils Thuerey + * + * This program is free software, distributed under the terms of the + * Apache License, Version 2.0 + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Accurate Viscous Free Surfaces for Buckling, Coiling, and Rotating Liquids + * Batty et al., SCA 2008 + * + ******************************************************************************/ + +#include "conjugategrad.h" +#include "general.h" +#include "grid.h" +#include "vectorbase.h" + +#include <chrono> + +#if OPENMP == 1 || TBB == 1 +# define ENABLE_PARALLEL 0 +#endif + +#if ENABLE_PARALLEL == 1 +# include <thread> +# include <algorithm> + +static const int manta_num_threads = std::thread::hardware_concurrency(); + +# define parallel_block \ + { \ + std::vector<std::thread> threads; \ + { + +# define do_parallel threads.push_back( std::thread([&]() { +# define do_end \ + } ) ); + +# define block_end \ + } \ + for (auto &thread : threads) { \ + thread.join(); \ + } \ + } + +#endif + +#define FOR_INT_IJK(num) \ + for (int k_off = 0; k_off < num; ++k_off) \ + for (int j_off = 0; j_off < num; ++j_off) \ + for (int i_off = 0; i_off < num; ++i_off) + +using namespace std; + +namespace Manta { + +//! Assumes phi0<0 and phi1>=0, phi2>=0, and phi3>=0 or vice versa. +//! In particular, phi0 must not equal any of phi1, phi2 or phi3. +static Real sortedTetFraction(Real phi0, Real phi1, Real phi2, Real phi3) +{ + return phi0 * phi0 * phi0 / ((phi0 - phi1) * (phi0 - phi2) * (phi0 - phi3)); +} + +//! Assumes phi0<0, phi1<0, and phi2>=0, and phi3>=0 or vice versa. +//! In particular, phi0 and phi1 must not equal any of phi2 and phi3. +static Real sortedPrismFraction(Real phi0, Real phi1, Real phi2, Real phi3) +{ + Real a = phi0 / (phi0 - phi2); + Real b = phi0 / (phi0 - phi3); + Real c = phi1 / (phi1 - phi3); + Real d = phi1 / (phi1 - phi2); + return a * b * (1 - d) + b * (1 - c) * d + c * d; +} + +Real volumeFraction(Real phi0, Real phi1, Real phi2, Real phi3) +{ + sort(phi0, phi1, phi2, phi3); + if (phi3 <= 0) + return 1; + else if (phi2 <= 0) + return 1 - sortedTetFraction(phi3, phi2, phi1, phi0); + else if (phi1 <= 0) + return sortedPrismFraction(phi0, phi1, phi2, phi3); + else if (phi0 <= 0) + return sortedTetFraction(phi0, phi1, phi2, phi3); + else + return 0; +} + +//! The average of the two possible decompositions of the cube into five tetrahedra. +Real volumeFraction(Real phi000, + Real phi100, + Real phi010, + Real phi110, + Real phi001, + Real phi101, + Real phi011, + Real phi111) +{ + return (volumeFraction(phi000, phi001, phi101, phi011) + + volumeFraction(phi000, phi101, phi100, phi110) + + volumeFraction(phi000, phi010, phi011, phi110) + + volumeFraction(phi101, phi011, phi111, phi110) + + 2 * volumeFraction(phi000, phi011, phi101, phi110) + + volumeFraction(phi100, phi101, phi001, phi111) + + volumeFraction(phi100, phi001, phi000, phi010) + + volumeFraction(phi100, phi110, phi111, phi010) + + volumeFraction(phi001, phi111, phi011, phi010) + + 2 * volumeFraction(phi100, phi111, phi001, phi010)) / + 12; +} + +//! Kernel loop over grid with 2x base resolution! + +struct KnEstimateVolumeFraction : public KernelBase { + KnEstimateVolumeFraction(Grid<Real> &volumes, + const Grid<Real> &phi, + const Vec3 &startCentre, + const Real dx) + : KernelBase(&volumes, 0), volumes(volumes), phi(phi), startCentre(startCentre), dx(dx) + { + runMessage(); + run(); + } + inline void op(int i, + int j, + int k, + Grid<Real> &volumes, + const Grid<Real> &phi, + const Vec3 &startCentre, + const Real dx) const + { + const Vec3 centre = startCentre + Vec3(i, j, k) * 0.5; + const Real offset = 0.5 * dx; + const int order = 2; + + Real phi000 = phi.getInterpolatedHi(centre + Vec3(-offset, -offset, -offset), order); + Real phi001 = phi.getInterpolatedHi(centre + Vec3(-offset, -offset, +offset), order); + Real phi010 = phi.getInterpolatedHi(centre + Vec3(-offset, +offset, -offset), order); + Real phi011 = phi.getInterpolatedHi(centre + Vec3(-offset, +offset, +offset), order); + Real phi100 = phi.getInterpolatedHi(centre + Vec3(+offset, -offset, -offset), order); + Real phi101 = phi.getInterpolatedHi(centre + Vec3(+offset, -offset, +offset), order); + Real phi110 = phi.getInterpolatedHi(centre + Vec3(+offset, +offset, -offset), order); + Real phi111 = phi.getInterpolatedHi(centre + Vec3(+offset, +offset, +offset), order); + + volumes(i, j, k) = volumeFraction( + phi000, phi100, phi010, phi110, phi001, phi101, phi011, phi111); + } + inline Grid<Real> &getArg0() + { + return volumes; + } + typedef Grid<Real> type0; + inline const Grid<Real> &getArg1() + { + return phi; + } + typedef Grid<Real> type1; + inline const Vec3 &getArg2() + { + return startCentre; + } + typedef Vec3 type2; + inline const Real &getArg3() + { + return dx; + } + typedef Real type3; + void runMessage() + { + debMsg("Executing kernel KnEstimateVolumeFraction ", 3); + debMsg("Kernel range" + << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ", + 4); + }; + void operator()(const tbb::blocked_range<IndexInt> &__r) const + { + const int _maxX = maxX; + const int _maxY = maxY; + if (maxZ > 1) { + for (int k = __r.begin(); k != (int)__r.end(); k++) + for (int j = 0; j < _maxY; j++) + for (int i = 0; i < _maxX; i++) + op(i, j, k, volumes, phi, startCentre, dx); + } + else { + const int k = 0; + for (int j = __r.begin(); j != (int)__r.end(); j++) + for (int i = 0; i < _maxX; i++) + op(i, j, k, volumes, phi, startCentre, dx); + } + } + void run() + { + if (maxZ > 1) + tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this); + else + tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this); + } + Grid<Real> &volumes; + const Grid<Real> φ + const Vec3 &startCentre; + const Real dx; +}; + +struct KnUpdateVolumeGrid : public KernelBase { + KnUpdateVolumeGrid(Grid<Real> &cVolLiquid, + Grid<Real> &uVolLiquid, + Grid<Real> &vVolLiquid, + Grid<Real> &wVolLiquid, + Grid<Real> &exVolLiquid, + Grid<Real> &eyVolLiquid, + Grid<Real> &ezVolLiquid, + const Grid<Real> &src) + : KernelBase(&cVolLiquid, 0), + cVolLiquid(cVolLiquid), + uVolLiquid(uVolLiquid), + vVolLiquid(vVolLiquid), + wVolLiquid(wVolLiquid), + exVolLiquid(exVolLiquid), + eyVolLiquid(eyVolLiquid), + ezVolLiquid(ezVolLiquid), + src(src) + { + runMessage(); + run(); + } + inline void op(int i, + int j, + int k, + Grid<Real> &cVolLiquid, + Grid<Real> &uVolLiquid, + Grid<Real> &vVolLiquid, + Grid<Real> &wVolLiquid, + Grid<Real> &exVolLiquid, + Grid<Real> &eyVolLiquid, + Grid<Real> &ezVolLiquid, + const Grid<Real> &src) const + { + // Work out c + cVolLiquid(i, j, k) = 0; + FOR_INT_IJK(2) + { + cVolLiquid(i, j, k) += src(2 * i + i_off, 2 * j + j_off, 2 * k + k_off); + } + cVolLiquid(i, j, k) /= 8; + + // Work out u + if (i >= 1) { + uVolLiquid(i, j, k) = 0; + int base_i = 2 * i - 1; + int base_j = 2 * j; + int base_k = 2 * k; + FOR_INT_IJK(2) + { + uVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off); + } + uVolLiquid(i, j, k) /= 8; + } + + // v + if (j >= 1) { + vVolLiquid(i, j, k) = 0; + int base_i = 2 * i; + int base_j = 2 * j - 1; + int base_k = 2 * k; + FOR_INT_IJK(2) + { + vVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off); + } + vVolLiquid(i, j, k) /= 8; + } + + // w + if (k >= 1) { + wVolLiquid(i, j, k) = 0; + int base_i = 2 * i; + int base_j = 2 * j; + int base_k = 2 * k - 1; + FOR_INT_IJK(2) + { + wVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off); + } + wVolLiquid(i, j, k) /= 8; + } + + // e-x + if (j >= 1 && k >= 1) { + exVolLiquid(i, j, k) = 0; + int base_i = 2 * i; + int base_j = 2 * j - 1; + int base_k = 2 * k - 1; + FOR_INT_IJK(2) + { + exVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off); + } + exVolLiquid(i, j, k) /= 8; + } + + // e-y + if (i >= 1 && k >= 1) { + eyVolLiquid(i, j, k) = 0; + int base_i = 2 * i - 1; + int base_j = 2 * j; + int base_k = 2 * k - 1; + FOR_INT_IJK(2) + { + eyVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off); + } + eyVolLiquid(i, j, k) /= 8; + } + + // e-z + if (i >= 1 && j >= 1) { + ezVolLiquid(i, j, k) = 0; + int base_i = 2 * i - 1; + int base_j = 2 * j - 1; + int base_k = 2 * k; + FOR_INT_IJK(2) + { + ezVolLiquid(i, j, k) += src(base_i + i_off, base_j + j_off, base_k + k_off); + } + ezVolLiquid(i, j, k) /= 8; + } + } + inline Grid<Real> &getArg0() + { + return cVolLiquid; + } + typedef Grid<Real> type0; + inline Grid<Real> &getArg1() + { + return uVolLiquid; + } + typedef Grid<Real> type1; + inline Grid<Real> &getArg2() + { + return vVolLiquid; + } + typedef Grid<Real> type2; + inline Grid<Real> &getArg3() + { + return wVolLiquid; + } + typedef Grid<Real> type3; + inline Grid<Real> &getArg4() + { + return exVolLiquid; + } + typedef Grid<Real> type4; + inline Grid<Real> &getArg5() + { + return eyVolLiquid; + } + typedef Grid<Real> type5; + inline Grid<Real> &getArg6() + { + return ezVolLiquid; + } + typedef Grid<Real> type6; + inline const Grid<Real> &getArg7() + { + return src; + } + typedef Grid<Real> type7; + void runMessage() + { + debMsg("Executing kernel KnUpdateVolumeGrid ", 3); + debMsg("Kernel range" + << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ", + 4); + }; + void operator()(const tbb::blocked_range<IndexInt> &__r) const + { + const int _maxX = maxX; + const int _maxY = maxY; + if (maxZ > 1) { + for (int k = __r.begin(); k != (int)__r.end(); k++) + for (int j = 0; j < _maxY; j++) + for (int i = 0; i < _maxX; i++) + op(i, + j, + k, + cVolLiquid, + uVolLiquid, + vVolLiquid, + wVolLiquid, + exVolLiquid, + eyVolLiquid, + ezVolLiquid, + src); + } + else { + const int k = 0; + for (int j = __r.begin(); j != (int)__r.end(); j++) + for (int i = 0; i < _maxX; i++) + op(i, + j, + k, + cVolLiquid, + uVolLiquid, + vVolLiquid, + wVolLiquid, + exVolLiquid, + eyVolLiquid, + ezVolLiquid, + src); + } + } + void run() + { + if (maxZ > 1) + tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this); + else + tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this); + } + Grid<Real> &cVolLiquid; + Grid<Real> &uVolLiquid; + Grid<Real> &vVolLiquid; + Grid<Real> &wVolLiquid; + Grid<Real> &exVolLiquid; + Grid<Real> &eyVolLiquid; + Grid<Real> &ezVolLiquid; + const Grid<Real> &src; +}; + +void computeWeights(const Grid<Real> &phi, + Grid<Real> &doubleSized, + Grid<Real> &cVolLiquid, + Grid<Real> &uVolLiquid, + Grid<Real> &vVolLiquid, + Grid<Real> &wVolLiquid, + Grid<Real> &exVolLiquid, + Grid<Real> &eyVolLiquid, + Grid<Real> &ezVolLiquid, + Real dx) +{ + KnEstimateVolumeFraction(doubleSized, phi, Vec3(0.25 * dx, 0.25 * dx, 0.25 * dx), 0.5 * dx); + KnUpdateVolumeGrid(cVolLiquid, + uVolLiquid, + vVolLiquid, + wVolLiquid, + exVolLiquid, + eyVolLiquid, + ezVolLiquid, + doubleSized); +} + +struct KnUpdateFaceStates : public KernelBase { + KnUpdateFaceStates(const FlagGrid &flags, + Grid<int> &uState, + Grid<int> &vState, + Grid<int> &wState) + : KernelBase(&flags, 0), flags(flags), uState(uState), vState(vState), wState(wState) + { + runMessage(); + run(); + } + inline void op(int i, + int j, + int k, + const FlagGrid &flags, + Grid<int> &uState, + Grid<int> &vState, + Grid<int> &wState) const + { + bool curObs = flags.isObstacle(i, j, k); + uState(i, j, k) = (i > 0 && !flags.isObstacle(i - 1, j, k) && !curObs) ? + FlagGrid::TypeFluid : + FlagGrid::TypeObstacle; + vState(i, j, k) = (j > 0 && !flags.isObstacle(i, j - 1, k) && !curObs) ? + FlagGrid::TypeFluid : + FlagGrid::TypeObstacle; + wState(i, j, k) = (k > 0 && !flags.isObstacle(i, j, k - 1) && !curObs) ? + FlagGrid::TypeFluid : + FlagGrid::TypeObstacle; + } + inline const FlagGrid &getArg0() + { + return flags; + } + typedef FlagGrid type0; + inline Grid<int> &getArg1() + { + return uState; + } + typedef Grid<int> type1; + inline Grid<int> &getArg2() + { + return vState; + } + typedef Grid<int> type2; + inline Grid<int> &getArg3() + { + return wState; + } + typedef Grid<int> type3; + void runMessage() + { + debMsg("Executing kernel KnUpdateFaceStates ", 3); + debMsg("Kernel range" + << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ", + 4); + }; + void operator()(const tbb::blocked_range<IndexInt> &__r) const + { + const int _maxX = maxX; + const int _maxY = maxY; + if (maxZ > 1) { + for (int k = __r.begin(); k != (int)__r.end(); k++) + for (int j = 0; j < _maxY; j++) + for (int i = 0; i < _maxX; i++) + op(i, j, k, flags, uState, vState, wState); + } + else { + const int k = 0; + for (int j = __r.begin(); j != (int)__r.end(); j++) + for (int i = 0; i < _maxX; i++) + op(i, j, k, flags, uState, vState, wState); + } + } + void run() + { + if (maxZ > 1) + tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this); + else + tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this); + } + const FlagGrid &flags; + Grid<int> &uState; + Grid<int> &vState; + Grid<int> &wState; +}; + +struct KnApplyVelocities : public KernelBase { + KnApplyVelocities(MACGrid &dst, + const Grid<int> &uState, + const Grid<int> &vState, + const Grid<int> &wState, + Grid<Real> &srcU, + Grid<Real> &srcV, + Grid<Real> &srcW) + : KernelBase(&dst, 0), + dst(dst), + uState(uState), + vState(vState), + wState(wState), + srcU(srcU), + srcV(srcV), + srcW(srcW) + { + runMessage(); + run(); + } + inline void op(int i, + int j, + int k, + MACGrid &dst, + const Grid<int> &uState, + const Grid<int> &vState, + const Grid<int> &wState, + Grid<Real> &srcU, + Grid<Real> &srcV, + Grid<Real> &srcW) const + { + dst(i, j, k).x = (uState(i, j, k) == FlagGrid::TypeFluid) ? srcU(i, j, k) : 0; + dst(i, j, k).y = (vState(i, j, k) == FlagGrid::TypeFluid) ? srcV(i, j, k) : 0; + if (dst.is3D()) + dst(i, j, k).z = (wState(i, j, k) == FlagGrid::TypeFluid) ? srcW(i, j, k) : 0; + } + inline MACGrid &getArg0() + { + return dst; + } + typedef MACGrid type0; + inline const Grid<int> &getArg1() + { + return uState; + } + typedef Grid<int> type1; + inline const Grid<int> &getArg2() + { + return vState; + } + typedef Grid<int> type2; + inline const Grid<int> &getArg3() + { + return wState; + } + typedef Grid<int> type3; + inline Grid<Real> &getArg4() + { + return srcU; + } + typedef Grid<Real> type4; + inline Grid<Real> &getArg5() + { + return srcV; + } + typedef Grid<Real> type5; + inline Grid<Real> &getArg6() + { + return srcW; + } + typedef Grid<Real> type6; + void runMessage() + { + debMsg("Executing kernel KnApplyVelocities ", 3); + debMsg("Kernel range" + << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ", + 4); + }; + void operator()(const tbb::blocked_range<IndexInt> &__r) const + { + const int _maxX = maxX; + const int _maxY = maxY; + if (maxZ > 1) { + for (int k = __r.begin(); k != (int)__r.end(); k++) + for (int j = 0; j < _maxY; j++) + for (int i = 0; i < _maxX; i++) + op(i, j, k, dst, uState, vState, wState, srcU, srcV, srcW); + } + else { + const int k = 0; + for (int j = __r.begin(); j != (int)__r.end(); j++) + for (int i = 0; i < _maxX; i++) + op(i, j, k, dst, uState, vState, wState, srcU, srcV, srcW); + } + } + void run() + { + if (maxZ > 1) + tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this); + else + tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this); + } + MACGrid &dst; + const Grid<int> &uState; + const Grid<int> &vState; + const Grid<int> &wState; + Grid<Real> &srcU; + Grid<Real> &srcV; + Grid<Real> &srcW; +}; + +void solveViscosity(const FlagGrid &flags, + MACGrid &vel, + Grid<Real> &cVolLiquid, + Grid<Real> &uVolLiquid, + Grid<Real> &vVolLiquid, + Grid<Real> &wVolLiquid, + Grid<Real> &exVolLiquid, + Grid<Real> &eyVolLiquid, + Grid<Real> &ezVolLiquid, + Grid<Real> &viscosity, + const Real dt, + const Real dx, + const Real cgAccuracy, + const Real cgMaxIterFac) +{ + const Real factor = dt * square(1.0 / dx); + const int maxIter = (int)(cgMaxIterFac * flags.getSize().max()) * (flags.is3D() ? 1 : 4); + GridCg<ApplyMatrixViscosityU> *uGcg; + GridCg<ApplyMatrixViscosityV> *vGcg; + GridCg<ApplyMatrixViscosityW> *wGcg; + + // Tmp grids for CG solve in U, V, W dimensions + FluidSolver *parent = flags.getParent(); + Grid<Real> uResidual(parent); + Grid<Real> vResidual(parent); + Grid<Real> wResidual(parent); + Grid<Real> uSearch(parent); + Grid<Real> vSearch(parent); + Grid<Real> wSearch(parent); + Grid<Real> uTmp(parent); + Grid<Real> vTmp(parent); + Grid<Real> wTmp(parent); + Grid<Real> uRhs(parent); + Grid<Real> vRhs(parent); + Grid<Real> wRhs(parent); + + // A matrix U grids + Grid<Real> uA0(parent); // diagonal elements in A matrix + Grid<Real> uAplusi(parent); // neighbor at i+1 + Grid<Real> uAplusj(parent); // neighbor at j+1 + Grid<Real> uAplusk(parent); // neighbor at k+1 + Grid<Real> uAminusi(parent); // neighbor at i-1 + Grid<Real> uAminusj(parent); // neighbor at j-1 + Grid<Real> uAminusk(parent); // neighbor at k-1 + Grid<Real> uAhelper1(parent); // additional helper grids for off diagonal elements + Grid<Real> uAhelper2(parent); + Grid<Real> uAhelper3(parent); + Grid<Real> uAhelper4(parent); + Grid<Real> uAhelper5(parent); + Grid<Real> uAhelper6(parent); + Grid<Real> uAhelper7(parent); + Grid<Real> uAhelper8(parent); + + // A matrix V grids + Grid<Real> vA0(parent); + Grid<Real> vAplusi(parent); + Grid<Real> vAplusj(parent); + Grid<Real> vAplusk(parent); + Grid<Real> vAminusi(parent); + Grid<Real> vAminusj(parent); + Grid<Real> vAminusk(parent); + Grid<Real> vAhelper1(parent); + Grid<Real> vAhelper2(parent); + Grid<Real> vAhelper3(parent); + Grid<Real> vAhelper4(parent); + Grid<Real> vAhelper5(parent); + Grid<Real> vAhelper6(parent); + Grid<Real> vAhelper7(parent); + Grid<Real> vAhelper8(parent); + + // A matrix W grids + Grid<Real> wA0(parent); + Grid<Real> wAplusi(parent); + Grid<Real> wAplusj(parent); + Grid<Real> wAplusk(parent); + Grid<Real> wAminusi(parent); + Grid<Real> wAminusj(parent); + Grid<Real> wAminusk(parent); + Grid<Real> wAhelper1(parent); + Grid<Real> wAhelper2(parent); + Grid<Real> wAhelper3(parent); + Grid<Real> wAhelper4(parent); + Grid<Real> wAhelper5(parent); + Grid<Real> wAhelper6(parent); + Grid<Real> wAhelper7(parent); + Grid<Real> wAhelper8(parent); + + // Solution grids for CG solvers + Grid<Real> uSolution(parent); + Grid<Real> vSolution(parent); + Grid<Real> wSolution(parent); + + // Save state of voxel face (fluid or obstacle) + Grid<int> uState(parent); + Grid<int> vState(parent); + Grid<int> wState(parent); + + // Save state of voxel face (fluid or obstacle) + KnUpdateFaceStates(flags, uState, vState, wState); + + // Shorter names for flags, we will use them often + int isFluid = FlagGrid::TypeFluid; + int isObstacle = FlagGrid::TypeObstacle; + + // Main viscosity loop: construct A matrices and rhs's in all dimensions + FOR_IJK_BND(flags, 1) + { + + // U-terms: 2u_xx+ v_xy +uyy + u_zz + w_xz + if (uState(i, j, k) == isFluid) { + + uRhs(i, j, k) = uVolLiquid(i, j, k) * vel(i, j, k).x; + uA0(i, j, k) = uVolLiquid(i, j, k); + + Real viscRight = viscosity(i, j, k); + Real viscLeft = viscosity(i - 1, j, k); + Real volRight = cVolLiquid(i, j, k); + Real volLeft = cVolLiquid(i - 1, j, k); + + Real viscTop = 0.25 * (viscosity(i - 1, j + 1, k) + viscosity(i - 1, j, k) + + viscosity(i, j + 1, k) + viscosity(i, j, k)); + Real viscBottom = 0.25 * (viscosity(i - 1, j, k) + viscosity(i - 1, j - 1, k) + + viscosity(i, j, k) + viscosity(i, j - 1, k)); + Real volTop = ezVolLiquid(i, j + 1, k); + Real volBottom = ezVolLiquid(i, j, k); + + Real viscFront = 0.25 * (viscosity(i - 1, j, k + 1) + viscosity(i - 1, j, k) + + viscosity(i, j, k + 1) + viscosity(i, j, k)); + Real viscBack = 0.25 * (viscosity(i - 1, j, k) + viscosity(i - 1, j, k - 1) + + viscosity(i, j, k) + viscosity(i, j, k - 1)); + Real volFront = eyVolLiquid(i, j, k + 1); + Real volBack = eyVolLiquid(i, j, k); + + Real factorRight = 2 * factor * viscRight * volRight; + Real factorLeft = 2 * factor * viscLeft * volLeft; + Real factorTop = factor * viscTop * volTop; + Real factorBottom = factor * viscBottom * volBottom; + Real factorFront = factor * viscFront * volFront; + Real factorBack = factor * viscBack * volBack; + + // u_x_right + uA0(i, j, k) += factorRight; + if (uState(i + 1, j, k) == isFluid) { + uAplusi(i, j, k) += -factorRight; + } + else if (uState(i + 1, j, k) == isObstacle) { + uRhs(i, j, k) -= -vel(i + 1, j, k).x * factorRight; + } + + // u_x_left + uA0(i, j, k) += factorLeft; + if (uState(i - 1, j, k) == isFluid) { + uAminusi(i, j, k) += -factorLeft; + } + else if (uState(i - 1, j, k) == isObstacle) { + uRhs(i, j, k) -= -vel(i - 1, j, k).x * factorLeft; + } + + // u_y_top + uA0(i, j, k) += factorTop; + if (uState(i, j + 1, k) == isFluid) { + uAplusj(i, j, k) += -factorTop; + } + else if (uState(i, j + 1, k) == isObstacle) { + uRhs(i, j, k) -= -vel(i, j + 1, k).x * factorTop; + } + + // u_y_bottom + uA0(i, j, k) += factorBottom; + if (uState(i, j - 1, k) == isFluid) { + uAminusj(i, j, k) += -factorBottom; + } + else if (uState(i, j - 1, k) == isObstacle) { + uRhs(i, j, k) -= -vel(i, j - 1, k).x * factorBottom; + } + + // u_z_front + uA0(i, j, k) += factorFront; + if (uState(i, j, k + 1) == isFluid) { + uAplusk(i, j, k) += -factorFront; + } + else if (uState(i, j, k + 1) == isObstacle) { + uRhs(i, j, k) -= -vel(i, j, k + 1).x * factorFront; + } + + // u_z_back + uA0(i, j, k) += factorBack; + if (uState(i, j, k - 1) == isFluid) { + uAminusk(i, j, k) += -factorBack; + } + else if (uState(i, j, k - 1) == isObstacle) { + uRhs(i, j, k) -= -vel(i, j, k - 1).x * factorBack; + } + + // v_x_top + if (vState(i, j + 1, k) == isFluid) { + uAhelper1(i, j, k) += -factorTop; + } + else if (vState(i, j + 1, k) == isObstacle) { + uRhs(i, j, k) -= -vel(i, j + 1, k).y * factorTop; + } + + if (vState(i - 1, j + 1, k) == isFluid) { + uAhelper2(i, j, k) += factorTop; + } + else if (vState(i - 1, j + 1, k) == isObstacle) { + uRhs(i, j, k) -= vel(i - 1, j + 1, k).y * factorTop; + } + + // v_x_bottom + if (vState(i, j, k) == isFluid) { + uAhelper3(i, j, k) += factorBottom; + } + else if (vState(i, j, k) == isObstacle) { + uRhs(i, j, k) -= vel(i, j, k).y * factorBottom; + } + + if (vState(i - 1, j, k) == isFluid) { + uAhelper4(i, j, k) += -factorBottom; + } + else if (vState(i - 1, j, k) == isObstacle) { + uRhs(i, j, k) -= -vel(i - 1, j, k).y * factorBottom; + } + + // w_x_front + if (wState(i, j, k + 1) == isFluid) { + uAhelper5(i, j, k) += -factorFront; + } + else if (wState(i, j, k + 1) == isObstacle) { + uRhs(i, j, k) -= -vel(i, j, k + 1).z * factorFront; + } + + if (wState(i - 1, j, k + 1) == isFluid) { + uAhelper6(i, j, k) += factorFront; + } + else if (wState(i - 1, j, k + 1) == isObstacle) { + uRhs(i, j, k) -= vel(i - 1, j, k + 1).z * factorFront; + } + + // w_x_back + if (wState(i, j, k) == isFluid) { + uAhelper7(i, j, k) += factorBack; + } + else if (wState(i, j, k) == isObstacle) { + uRhs(i, j, k) -= vel(i, j, k).z * factorBack; + } + + if (wState(i - 1, j, k) == isFluid) { + uAhelper8(i, j, k) += -factorBack; + } + else if (wState(i - 1, j, k) == isObstacle) { + uRhs(i, j, k) -= -vel(i - 1, j, k).z * factorBack; + } + } + + // V-terms: vxx + 2vyy + vzz + u_yx + w_yz + if (vState(i, j, k) == isFluid) { + + vRhs(i, j, k) = vVolLiquid(i, j, k) * vel(i, j, k).y; + vA0(i, j, k) = vVolLiquid(i, j, k); + + Real viscRight = 0.25 * (viscosity(i, j - 1, k) + viscosity(i + 1, j - 1, k) + + viscosity(i, j, k) + viscosity(i + 1, j, k)); + Real viscLeft = 0.25 * (viscosity(i, j - 1, k) + viscosity(i - 1, j - 1, k) + + viscosity(i, j, k) + viscosity(i - 1, j, k)); + Real volRight = ezVolLiquid(i + 1, j, k); + Real volLeft = ezVolLiquid(i, j, k); + + Real viscTop = viscosity(i, j, k); + Real viscBottom = viscosity(i, j - 1, k); + Real volTop = cVolLiquid(i, j, k); + Real volBottom = cVolLiquid(i, j - 1, k); + + Real viscFront = 0.25 * (viscosity(i, j - 1, k) + viscosity(i, j - 1, k + 1) + + viscosity(i, j, k) + viscosity(i, j, k + 1)); + Real viscBack = 0.25 * (viscosity(i, j - 1, k) + viscosity(i, j - 1, k - 1) + + viscosity(i, j, k) + viscosity(i, j, k - 1)); + Real volFront = exVolLiquid(i, j, k + 1); + Real volBack = exVolLiquid(i, j, k); + + Real factorRight = factor * viscRight * volRight; + Real factorLeft = factor * viscLeft * volLeft; + Real factorTop = 2 * factor * viscTop * volTop; + Real factorBottom = 2 * factor * viscBottom * volBottom; + Real factorFront = factor * viscFront * volFront; + Real factorBack = factor * viscBack * volBack; + + // v_x_right + vA0(i, j, k) += factorRight; + if (vState(i + 1, j, k) == isFluid) { + vAplusi(i, j, k) += -factorRight; + } + else if (vState(i + 1, j, k) == isObstacle) { + vRhs(i, j, k) -= -vel(i + 1, j, k).y * factorRight; + } + + // v_x_left + vA0(i, j, k) += factorLeft; + if (vState(i - 1, j, k) == isFluid) { + vAminusi(i, j, k) += -factorLeft; + } + else if (vState(i - 1, j, k) == isObstacle) { + vRhs(i, j, k) -= -vel(i - 1, j, k).y * factorLeft; + } + + // vy_top + vA0(i, j, k) += factorTop; + if (vState(i, j + 1, k) == isFluid) { + vAplusj(i, j, k) += -factorTop; + } + else if (vState(i, j + 1, k) == isObstacle) { + vRhs(i, j, k) -= -vel(i, j + 1, k).y * factorTop; + } + + // vy_bottom + vA0(i, j, k) += factorBottom; + if (vState(i, j - 1, k) == isFluid) { + vAminusj(i, j, k) += -factorBottom; + } + else if (vState(i, j - 1, k) == isObstacle) { + vRhs(i, j, k) -= -vel(i, j - 1, k).y * factorBottom; + } + + // v_z_front + vA0(i, j, k) += factorFront; + if (vState(i, j, k + 1) == isFluid) { + vAplusk(i, j, k) += -factorFront; + } + else if (vState(i, j, k + 1) == isObstacle) { + vRhs(i, j, k) -= -vel(i, j, k + 1).y * factorFront; + } + + // v_z_back + vA0(i, j, k) += factorBack; + if (vState(i, j, k - 1) == isFluid) { + vAminusk(i, j, k) += -factorBack; + } + else if (vState(i, j, k - 1) == isObstacle) { + vRhs(i, j, k) -= -vel(i, j, k - 1).y * factorBack; + } + + // u_y_right + if (uState(i + 1, j, k) == isFluid) { + vAhelper1(i, j, k) += -factorRight; + } + else if (uState(i + 1, j, k) == isObstacle) { + vRhs(i, j, k) -= -vel(i + 1, j, k).x * factorRight; + } + + if (uState(i + 1, j - 1, k) == isFluid) { + vAhelper2(i, j, k) += factorRight; + } + else if (uState(i + 1, j - 1, k) == isObstacle) { + vRhs(i, j, k) -= vel(i + 1, j - 1, k).x * factorRight; + } + + // u_y_left + if (uState(i, j, k) == isFluid) { + vAhelper3(i, j, k) += factorLeft; + } + else if (uState(i, j, k) == isObstacle) { + vRhs(i, j, k) -= vel(i, j, k).x * factorLeft; + } + + if (uState(i, j - 1, k) == isFluid) { + vAhelper4(i, j, k) += -factorLeft; + } + else if (uState(i, j - 1, k) == isObstacle) { + vRhs(i, j, k) -= -vel(i, j - 1, k).x * factorLeft; + } + + // w_y_front + if (wState(i, j, k + 1) == isFluid) { + vAhelper5(i, j, k) += -factorFront; + } + else if (wState(i, j, k + 1) == isObstacle) { + vRhs(i, j, k) -= -vel(i, j, k + 1).z * factorFront; + } + + if (wState(i, j - 1, k + 1) == isFluid) { + vAhelper6(i, j, k) += factorFront; + } + else if (wState(i, j - 1, k + 1) == isObstacle) { + vRhs(i, j, k) -= vel(i, j - 1, k + 1).z * factorFront; + } + + // w_y_back + if (wState(i, j, k) == isFluid) { + vAhelper7(i, j, k) += factorBack; + } + else if (wState(i, j, k) == isObstacle) { + vRhs(i, j, k) -= vel(i, j, k).z * factorBack; + } + + if (wState(i, j - 1, k) == isFluid) { + vAhelper8(i, j, k) += -factorBack; + } + else if (wState(i, j - 1, k) == isObstacle) { + vRhs(i, j, k) -= -vel(i, j - 1, k).z * factorBack; + } + } + + // W-terms: wxx+ wyy+ 2wzz + u_zx + v_zy + if (wState(i, j, k) == isFluid) { + + wRhs(i, j, k) = wVolLiquid(i, j, k) * vel(i, j, k).z; + wA0(i, j, k) = wVolLiquid(i, j, k); + + Real viscRight = 0.25 * (viscosity(i, j, k) + viscosity(i, j, k - 1) + + viscosity(i + 1, j, k) + viscosity(i + 1, j, k - 1)); + Real viscLeft = 0.25 * (viscosity(i, j, k) + viscosity(i, j, k - 1) + + viscosity(i - 1, j, k) + viscosity(i - 1, j, k - 1)); + Real volRight = eyVolLiquid(i + 1, j, k); + Real volLeft = eyVolLiquid(i, j, k); + + Real viscTop = 0.25 * (viscosity(i, j, k) + viscosity(i, j, k - 1) + viscosity(i, j + 1, k) + + viscosity(i, j + 1, k - 1)); + ; + Real viscBottom = 0.25 * (viscosity(i, j, k) + viscosity(i, j, k - 1) + + viscosity(i, j - 1, k) + viscosity(i, j - 1, k - 1)); + ; + Real volTop = exVolLiquid(i, j + 1, k); + Real volBottom = exVolLiquid(i, j, k); + + Real viscFront = viscosity(i, j, k); + Real viscBack = viscosity(i, j, k - 1); + Real volFront = cVolLiquid(i, j, k); + Real volBack = cVolLiquid(i, j, k - 1); + + Real factorRight = factor * viscRight * volRight; + Real factorLeft = factor * viscLeft * volLeft; + Real factorTop = factor * viscTop * volTop; + Real factorBottom = factor * viscBottom * volBottom; + Real factorFront = 2 * factor * viscFront * volFront; + Real factorBack = 2 * factor * viscBack * volBack; + + // w_x_right + wA0(i, j, k) += factorRight; + if (wState(i + 1, j, k) == isFluid) { + wAplusi(i, j, k) += -factorRight; + } + else if (wState(i + 1, j, k) == isObstacle) { + wRhs(i, j, k) -= -vel(i + 1, j, k).z * factorRight; + } + + // w_x_left + wA0(i, j, k) += factorLeft; + if (wState(i - 1, j, k) == isFluid) { + wAminusi(i, j, k) += -factorLeft; + } + else if (wState(i - 1, j, k) == isObstacle) { + wRhs(i, j, k) -= -vel(i - 1, j, k).z * factorLeft; + } + + // w_y_top + wA0(i, j, k) += factorTop; + if (wState(i, j + 1, k) == isFluid) { + wAplusj(i, j, k) += -factorTop; + } + else if (wState(i, j + 1, k) == isObstacle) { + wRhs(i, j, k) -= -vel(i, j + 1, k).z * factorTop; + } + + // w_y_bottom + wA0(i, j, k) += factorBottom; + if (wState(i, j - 1, k) == isFluid) { + wAminusj(i, j, k) += -factorBottom; + } + else if (wState(i, j - 1, k) == isObstacle) { + wRhs(i, j, k) -= -vel(i, j - 1, k).z * factorBottom; + } + + // w_z_front + wA0(i, j, k) += factorFront; + if (wState(i, j, k + 1) == isFluid) { + wAplusk(i, j, k) += -factorFront; + } + else if (wState(i, j, k + 1) == isObstacle) { + wRhs(i, j, k) -= -vel(i, j, k + 1).z * factorFront; + } + + // w_z_back + wA0(i, j, k) += factorBack; + if (wState(i, j, k - 1) == isFluid) { + wAminusk(i, j, k) += -factorBack; + } + else if (wState(i, j, k - 1) == isObstacle) { + wRhs(i, j, k) -= -vel(i, j, k - 1).z * factorBack; + } + + // u_z_right + if (uState(i + 1, j, k) == isFluid) { + wAhelper1(i, j, k) += -factorRight; + } + else if (uState(i + 1, j, k) == isObstacle) { + wRhs(i, j, k) -= -vel(i + 1, j, k).x * factorRight; + } + + if (uState(i + 1, j, k - 1) == isFluid) { + wAhelper2(i, j, k) += factorRight; + } + else if (uState(i + 1, j, k - 1) == isObstacle) { + wRhs(i, j, k) -= vel(i + 1, j, k - 1).x * factorRight; + } + + // u_z_left + if (uState(i, j, k) == isFluid) { + wAhelper3(i, j, k) += factorLeft; + } + else if (uState(i, j, k) == isObstacle) { + wRhs(i, j, k) -= vel(i, j, k).x * factorLeft; + } + + if (uState(i, j, k - 1) == isFluid) { + wAhelper4(i, j, k) += -factorLeft; + } + else if (uState(i, j, k - 1) == isObstacle) { + wRhs(i, j, k) -= -vel(i, j, k - 1).x * factorLeft; + } + + // v_z_top + if (vState(i, j + 1, k) == isFluid) { + wAhelper5(i, j, k) += -factorTop; + } + else if (vState(i, j + 1, k) == isObstacle) { + wRhs(i, j, k) -= -vel(i, j + 1, k).y * factorTop; + } + + if (vState(i, j + 1, k - 1) == isFluid) { + wAhelper6(i, j, k) += factorTop; + } + else if (vState(i, j + 1, k - 1) == isObstacle) { + wRhs(i, j, k) -= vel(i, j + 1, k - 1).y * factorTop; + } + + // v_z_bottom + if (vState(i, j, k) == isFluid) { + wAhelper7(i, j, k) += factorBottom; + } + else if (vState(i, j, k) == isObstacle) { + wRhs(i, j, k) -= vel(i, j, k).y * factorBottom; + } + + if (vState(i, j, k - 1) == isFluid) { + wAhelper8(i, j, k) += -factorBottom; + } + else if (vState(i, j, k - 1) == isObstacle) { + wRhs(i, j, k) -= -vel(i, j, k - 1).y * factorBottom; + } + } + } + + // CG solver for U + if (flags.is3D()) { + vector<Grid<Real> *> uMatA{&uA0, + &uAplusi, + &uAplusj, + &uAplusk, + &uAminusi, + &uAminusj, + &uAminusk, + &uAhelper1, + &uAhelper2, + &uAhelper3, + &uAhelper4, + &uAhelper5, + &uAhelper6, + &uAhelper7, + &uAhelper8}; + vector<Grid<Real> *> uVecRhs{&vRhs, &wRhs}; + uGcg = new GridCg<ApplyMatrixViscosityU>( + uSolution, uRhs, uResidual, uSearch, flags, uTmp, uMatA, uVecRhs); + } + else { + errMsg("2D Matrix application not yet supported in viscosity solver"); + } + + // CG solver for V + if (flags.is3D()) { + vector<Grid<Real> *> vMatA{&vA0, + &vAplusi, + &vAplusj, + &vAplusk, + &vAminusi, + &vAminusj, + &vAminusk, + &vAhelper1, + &vAhelper2, + &vAhelper3, + &vAhelper4, + &vAhelper5, + &vAhelper6, + &vAhelper7, + &vAhelper8}; + vector<Grid<Real> *> vVecRhs{&uRhs, &wRhs}; + vGcg = new GridCg<ApplyMatrixViscosityV>( + vSolution, vRhs, vResidual, vSearch, flags, vTmp, vMatA, vVecRhs); + } + else { + errMsg("2D Matrix application not yet supported in viscosity solver"); + } + + // CG solver for W + if (flags.is3D()) { + vector<Grid<Real> *> wMatA{&wA0, + &wAplusi, + &wAplusj, + &wAplusk, + &wAminusi, + &wAminusj, + &wAminusk, + &wAhelper1, + &wAhelper2, + &wAhelper3, + &wAhelper4, + &wAhelper5, + &wAhelper6, + &wAhelper7, + &wAhelper8}; + vector<Grid<Real> *> wVecRhs{&uRhs, &vRhs}; + wGcg = new GridCg<ApplyMatrixViscosityW>( + wSolution, wRhs, wResidual, wSearch, flags, wTmp, wMatA, wVecRhs); + } + else { + errMsg("2D Matrix application not yet supported in viscosity solver"); + } + + // Same accuracy for all dimensions + uGcg->setAccuracy(cgAccuracy); + vGcg->setAccuracy(cgAccuracy); + wGcg->setAccuracy(cgAccuracy); + + // CG solve. Preconditioning not supported yet. Instead, U, V, W can optionally be solved in + // parallel. + for (int uIter = 0, vIter = 0, wIter = 0; uIter < maxIter || vIter < maxIter || wIter < maxIter; + uIter++, vIter++, wIter++) { +#if ENABLE_PARALLEL == 1 + parallel_block do_parallel +#endif + if (uIter < maxIter && !uGcg->iterate()) uIter = maxIter; +#if ENABLE_PARALLEL == 1 + do_end do_parallel +#endif + if (vIter < maxIter && !vGcg->iterate()) vIter = maxIter; +#if ENABLE_PARALLEL == 1 + do_end do_parallel +#endif + if (wIter < maxIter && !wGcg->iterate()) wIter = maxIter; +#if ENABLE_PARALLEL == 1 + do_end block_end +#endif + + // Make sure that next CG iteration has updated rhs grids + uRhs.copyFrom(uSearch); + vRhs.copyFrom(vSearch); + wRhs.copyFrom(wSearch); + } + debMsg( + "Viscosity::solveViscosity done. " + "Iterations (u,v,w): (" + << uGcg->getIterations() << "," << vGcg->getIterations() << "," << wGcg->getIterations() + << "), " + "Residual (u,v,w): (" + << uGcg->getResNorm() << "," << vGcg->getResNorm() << "," << wGcg->getResNorm() << ")", + 2); + + delete uGcg; + delete vGcg; + delete wGcg; + + // Apply solutions to global velocity grid + KnApplyVelocities(vel, uState, vState, wState, uSolution, vSolution, wSolution); +} + +//! To use the viscosity plugin, scenes must call this function before solving pressure. +//! Note that the 'volumes' grid uses 2x the base resolution + +void applyViscosity(const FlagGrid &flags, + const Grid<Real> &phi, + MACGrid &vel, + Grid<Real> &volumes, + Grid<Real> &viscosity, + const Real cgAccuracy = 1e-9, + const Real cgMaxIterFac = 1.5) +{ + const Real dx = flags.getParent()->getDx(); + const Real dt = flags.getParent()->getDt(); + + // Reserve temp grids for volume weight calculation + FluidSolver *parent = flags.getParent(); + Grid<Real> cVolLiquid(parent); + Grid<Real> uVolLiquid(parent); + Grid<Real> vVolLiquid(parent); + Grid<Real> wVolLiquid(parent); + Grid<Real> exVolLiquid(parent); + Grid<Real> eyVolLiquid(parent); + Grid<Real> ezVolLiquid(parent); + + // Ensure final weight grid gets cleared at every step + volumes.clear(); + + // Save viscous fluid volume in double-sized volumes grid + computeWeights(phi, + volumes, + cVolLiquid, + uVolLiquid, + vVolLiquid, + wVolLiquid, + exVolLiquid, + eyVolLiquid, + ezVolLiquid, + dx); + + // Set up A matrix and rhs. Solve with CG. Update velocity grid. + solveViscosity(flags, + vel, + cVolLiquid, + uVolLiquid, + vVolLiquid, + wVolLiquid, + exVolLiquid, + eyVolLiquid, + ezVolLiquid, + viscosity, + dt, + dx, + cgAccuracy, + cgMaxIterFac); +} +static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds) +{ + try { + PbArgs _args(_linargs, _kwds); + FluidSolver *parent = _args.obtainParent(); + bool noTiming = _args.getOpt<bool>("notiming", -1, 0); + pbPreparePlugin(parent, "applyViscosity", !noTiming); + PyObject *_retval = nullptr; + { + ArgLocker _lock; + const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock); + const Grid<Real> &phi = *_args.getPtr<Grid<Real>>("phi", 1, &_lock); + MACGrid &vel = *_args.getPtr<MACGrid>("vel", 2, &_lock); + Grid<Real> &volumes = *_args.getPtr<Grid<Real>>("volumes", 3, &_lock); + Grid<Real> &viscosity = *_args.getPtr<Grid<Real>>("viscosity", 4, &_lock); + const Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 5, 1e-9, &_lock); + const Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 6, 1.5, &_lock); + _retval = getPyNone(); + applyViscosity(flags, phi, vel, volumes, viscosity, cgAccuracy, cgMaxIterFac); + _args.check(); + } + pbFinalizePlugin(parent, "applyViscosity", !noTiming); + return _retval; + } + catch (std::exception &e) { + pbSetError("applyViscosity", e.what()); + return 0; + } +} +static const Pb::Register _RP_applyViscosity("", "applyViscosity", _W_0); +extern "C" { +void PbRegister_applyViscosity() +{ + KEEP_UNUSED(_RP_applyViscosity); +} +} + +} // namespace Manta + +#if ENABLE_PARALLEL == 1 + +# undef parallel_block +# undef do_parallel +# undef do_end +# undef block_end +# undef parallel_for +# undef parallel_end + +#endif diff --git a/extern/mantaflow/preprocessed/plugin/vortexplugins.cpp b/extern/mantaflow/preprocessed/plugin/vortexplugins.cpp index d5d7d597a7c..18222c4ccda 100644 --- a/extern/mantaflow/preprocessed/plugin/vortexplugins.cpp +++ b/extern/mantaflow/preprocessed/plugin/vortexplugins.cpp @@ -576,8 +576,10 @@ void VICintegration(VortexSheetMesh &mesh, // prepare CG solver const int maxIter = (int)(cgMaxIterFac * vel.getSize().max()); + vector<Grid<Real> *> matA{&A0, &Ai, &Aj, &Ak}; + GridCgInterface *gcg = new GridCg<ApplyMatrix>( - solution, rhs, residual, search, flags, temp1, &A0, &Ai, &Aj, &Ak); + solution, rhs, residual, search, flags, temp1, matA); gcg->setAccuracy(cgAccuracy); gcg->setUseL2Norm(true); gcg->setICPreconditioner( diff --git a/extern/mantaflow/preprocessed/plugin/waves.cpp b/extern/mantaflow/preprocessed/plugin/waves.cpp index 6fb3b16c742..53c56b8c506 100644 --- a/extern/mantaflow/preprocessed/plugin/waves.cpp +++ b/extern/mantaflow/preprocessed/plugin/waves.cpp @@ -423,10 +423,15 @@ void cgSolveWE(const FlagGrid &flags, const int maxIter = (int)(cgMaxIterFac * flags.getSize().max()) * (flags.is3D() ? 1 : 4); GridCgInterface *gcg; - if (flags.is3D()) - gcg = new GridCg<ApplyMatrix>(out, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak); - else - gcg = new GridCg<ApplyMatrix2D>(out, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak); + vector<Grid<Real> *> matA{&A0, &Ai, &Aj}; + + if (flags.is3D()) { + matA.push_back(&Ak); + gcg = new GridCg<ApplyMatrix>(out, rhs, residual, search, flags, tmp, matA); + } + else { + gcg = new GridCg<ApplyMatrix2D>(out, rhs, residual, search, flags, tmp, matA); + } gcg->setAccuracy(cgAccuracy); diff --git a/extern/mantaflow/preprocessed/registration.cpp b/extern/mantaflow/preprocessed/registration.cpp index d5dae479f0e..fd32463b95f 100644 --- a/extern/mantaflow/preprocessed/registration.cpp +++ b/extern/mantaflow/preprocessed/registration.cpp @@ -145,6 +145,7 @@ extern void PbRegister_flipComputeSurfaceNormals(); extern void PbRegister_flipUpdateNeighborRatio(); extern void PbRegister_particleSurfaceTurbulence(); extern void PbRegister_debugCheckParts(); +extern void PbRegister_applyViscosity(); extern void PbRegister_markAsFixed(); extern void PbRegister_texcoordInflow(); extern void PbRegister_meshSmokeInflow(); @@ -342,6 +343,7 @@ void MantaEnsureRegistration() PbRegister_flipUpdateNeighborRatio(); PbRegister_particleSurfaceTurbulence(); PbRegister_debugCheckParts(); + PbRegister_applyViscosity(); PbRegister_markAsFixed(); PbRegister_texcoordInflow(); PbRegister_meshSmokeInflow(); |