diff options
24 files changed, 2148 insertions, 172 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(); diff --git a/intern/mantaflow/intern/MANTA_main.cpp b/intern/mantaflow/intern/MANTA_main.cpp index d49915e2452..d5d41e71f22 100644 --- a/intern/mantaflow/intern/MANTA_main.cpp +++ b/intern/mantaflow/intern/MANTA_main.cpp @@ -72,6 +72,7 @@ MANTA::MANTA(int *res, FluidModifierData *fmd) mUsingFractions = (fds->flags & FLUID_DOMAIN_USE_FRACTIONS) && mUsingLiquid; mUsingMesh = (fds->flags & FLUID_DOMAIN_USE_MESH) && mUsingLiquid; mUsingDiffusion = (fds->flags & FLUID_DOMAIN_USE_DIFFUSION) && mUsingLiquid; + mUsingViscosity = (fds->flags & FLUID_DOMAIN_USE_VISCOSITY) && mUsingLiquid; mUsingMVel = (fds->flags & FLUID_DOMAIN_USE_SPEED_VECTORS) && mUsingLiquid; mUsingGuiding = (fds->flags & FLUID_DOMAIN_USE_GUIDE); mUsingDrops = (fds->particle_type & FLUID_DOMAIN_PARTICLE_SPRAY) && mUsingLiquid; @@ -221,6 +222,10 @@ MANTA::MANTA(int *res, FluidModifierData *fmd) initSuccess &= initLiquidMesh(); } + if (mUsingViscosity) { + initSuccess &= initLiquidViscosity(); + } + if (mUsingDiffusion) { initSuccess &= initCurvature(); } @@ -440,6 +445,17 @@ bool MANTA::initLiquidMesh(FluidModifierData *fmd) return runPythonString(pythonCommands); } +bool MANTA::initLiquidViscosity(FluidModifierData *fmd) +{ + vector<string> pythonCommands; + string tmpString = fluid_variables_viscosity + fluid_solver_viscosity + liquid_alloc_viscosity; + string finalString = parseScript(tmpString, fmd); + pythonCommands.push_back(finalString); + + mUsingViscosity = true; + return runPythonString(pythonCommands); +} + bool MANTA::initCurvature(FluidModifierData *fmd) { std::vector<std::string> pythonCommands; @@ -871,8 +887,10 @@ void MANTA::initializeRNAMap(FluidModifierData *fmd) mRNAMap["CACHE_DIR"] = cacheDirectory; mRNAMap["COMPRESSION_OPENVDB"] = vdbCompressionMethod; mRNAMap["PRECISION_OPENVDB"] = vdbPrecisionHalf; - mRNAMap["CLIP_OPENVDB"] = to_string(fds->clipping); + mRNAMap["CLIP_OPENVDB"] = to_string(fds->clipping); mRNAMap["PP_PARTICLE_MAXIMUM"] = to_string(fds->sys_particle_maximum); + mRNAMap["USING_VISCOSITY"] = getBooleanString(fds->flags & FLUID_DOMAIN_USE_VISCOSITY); + mRNAMap["VISCOSITY_VALUE"] = to_string(fds->viscosity_value); /* Fluid object names. */ mRNAMap["NAME_FLAGS"] = FLUID_NAME_FLAGS; @@ -1728,6 +1746,7 @@ bool MANTA::exportLiquidScript(FluidModifierData *fmd) bool guiding = fds->active_fields & FLUID_DOMAIN_ACTIVE_GUIDE; bool invel = fds->active_fields & FLUID_DOMAIN_ACTIVE_INVEL; bool outflow = fds->active_fields & FLUID_DOMAIN_ACTIVE_OUTFLOW; + bool viscosity = fds->flags & FLUID_DOMAIN_USE_VISCOSITY; string manta_script; @@ -1742,6 +1761,8 @@ bool MANTA::exportLiquidScript(FluidModifierData *fmd) manta_script += fluid_variables_particles + liquid_variables_particles; if (guiding) manta_script += fluid_variables_guiding; + if (viscosity) + manta_script += fluid_variables_viscosity; /* Solvers. */ manta_script += header_solvers + fluid_solver; @@ -1751,6 +1772,8 @@ bool MANTA::exportLiquidScript(FluidModifierData *fmd) manta_script += fluid_solver_particles; if (guiding) manta_script += fluid_solver_guiding; + if (viscosity) + manta_script += fluid_solver_viscosity; /* Grids. */ manta_script += header_grids + fluid_alloc + liquid_alloc; @@ -1768,6 +1791,8 @@ bool MANTA::exportLiquidScript(FluidModifierData *fmd) manta_script += fluid_alloc_invel; if (outflow) manta_script += fluid_alloc_outflow; + if (viscosity) + manta_script += liquid_alloc_viscosity; /* Domain init. */ manta_script += header_gridinit + liquid_init_phi; diff --git a/intern/mantaflow/intern/MANTA_main.h b/intern/mantaflow/intern/MANTA_main.h index 163b168e43d..7129a545243 100644 --- a/intern/mantaflow/intern/MANTA_main.h +++ b/intern/mantaflow/intern/MANTA_main.h @@ -67,6 +67,7 @@ struct MANTA { bool initColorsHigh(struct FluidModifierData *fmd = nullptr); bool initLiquid(FluidModifierData *fmd = nullptr); bool initLiquidMesh(FluidModifierData *fmd = nullptr); + bool initLiquidViscosity(FluidModifierData *fmd = nullptr); bool initObstacle(FluidModifierData *fmd = nullptr); bool initCurvature(FluidModifierData *fmd = nullptr); bool initGuiding(FluidModifierData *fmd = nullptr); @@ -757,6 +758,7 @@ struct MANTA { bool mUsingNoise; bool mUsingMesh; bool mUsingDiffusion; + bool mUsingViscosity; bool mUsingMVel; bool mUsingLiquid; bool mUsingSmoke; diff --git a/intern/mantaflow/intern/strings/fluid_script.h b/intern/mantaflow/intern/strings/fluid_script.h index 8c9025dd435..73e8552d260 100644 --- a/intern/mantaflow/intern/strings/fluid_script.h +++ b/intern/mantaflow/intern/strings/fluid_script.h @@ -79,6 +79,11 @@ const std::string fluid_solver_guiding = mantaMsg('Solver guiding')\n\ sg$ID$ = Solver(name='solver_guiding$ID$', gridSize=gs_sg$ID$)\n"; +const std::string fluid_solver_viscosity = + "\n\ +mantaMsg('Solver viscosity')\n\ +sv$ID$ = Solver(name='solver_viscosity$ID$', gridSize=gs_sv$ID$, dim=dim_s$ID$)\n"; + ////////////////////////////////////////////////////////////////////// // VARIABLES ////////////////////////////////////////////////////////////////////// @@ -133,7 +138,7 @@ end_frame_s$ID$ = $END_FRAME$\n\ \n\ # Fluid diffusion / viscosity\n\ domainSize_s$ID$ = $FLUID_DOMAIN_SIZE$ # longest domain side in meters\n\ -viscosity_s$ID$ = $FLUID_VISCOSITY$ / (domainSize_s$ID$*domainSize_s$ID$) # kinematic viscosity in m^2/s\n\ +kinViscosity_s$ID$ = $FLUID_VISCOSITY$ / (domainSize_s$ID$*domainSize_s$ID$) # kinematic viscosity in m^2/s\n\ \n\ # Factors to convert Blender units to Manta units\n\ ratioMetersToRes_s$ID$ = float(domainSize_s$ID$) / float(res_s$ID$) # [meters / cells]\n\ @@ -199,6 +204,10 @@ tau_sg$ID$ = 1.0\n\ sigma_sg$ID$ = 0.99/tau_sg$ID$\n\ theta_sg$ID$ = 1.0\n"; +const std::string fluid_variables_viscosity = + "\n\ +gs_sv$ID$ = vec3($RESX$*2, $RESY$*2, $RESZ$*2)\n"; + const std::string fluid_with_obstacle = "\n\ using_obstacle_s$ID$ = True\n"; diff --git a/intern/mantaflow/intern/strings/liquid_script.h b/intern/mantaflow/intern/strings/liquid_script.h index 8e0d113d680..c44727bd47e 100644 --- a/intern/mantaflow/intern/strings/liquid_script.h +++ b/intern/mantaflow/intern/strings/liquid_script.h @@ -40,6 +40,8 @@ radiusFactor_s$ID$ = $PARTICLE_RADIUS$\n\ using_mesh_s$ID$ = $USING_MESH$\n\ using_final_mesh_s$ID$ = $USING_IMPROVED_MESH$\n\ using_fractions_s$ID$ = $USING_FRACTIONS$\n\ +using_apic_s$ID$ = $USING_APIC$\n\ +using_viscosity_s$ID$ = $USING_VISCOSITY$\n\ fracThreshold_s$ID$ = $FRACTIONS_THRESHOLD$\n\ fracDistance_s$ID$ = $FRACTIONS_DISTANCE$\n\ flipRatio_s$ID$ = $FLIP_RATIO$\n\ @@ -51,7 +53,7 @@ smoothenNeg_s$ID$ = $MESH_SMOOTHEN_NEG$\n\ randomness_s$ID$ = $PARTICLE_RANDOMNESS$\n\ surfaceTension_s$ID$ = $LIQUID_SURFACE_TENSION$\n\ maxSysParticles_s$ID$ = $PP_PARTICLE_MAXIMUM$\n\ -using_apic_s$ID$ = $USING_APIC$\n"; +viscosityValue_s$ID$ = $VISCOSITY_VALUE$\n"; const std::string liquid_variables_particles = "\n\ @@ -135,6 +137,13 @@ liquid_mesh_dict_s$ID$ = { 'lMesh' : mesh_sm$ID$ }\n\ if using_speedvectors_s$ID$:\n\ liquid_meshvel_dict_s$ID$ = { 'lVelMesh' : mVel_mesh$ID$ }\n"; +const std::string liquid_alloc_viscosity = + "\n\ +# Viscosity grids\n\ +volumes_s$ID$ = sv$ID$.create(RealGrid)\n\ +viscosity_s$ID$ = s$ID$.create(RealGrid)\n\ +viscosity_s$ID$.setConst(viscosityValue_s$ID$)\n"; + const std::string liquid_alloc_curvature = "\n\ mantaMsg('Liquid alloc curvature')\n\ @@ -306,7 +315,7 @@ def liquid_step_$ID$():\n\ if using_diffusion_s$ID$:\n\ mantaMsg('Viscosity')\n\ # diffusion param for solve = const * dt / dx^2\n\ - alphaV = viscosity_s$ID$ * s$ID$.timestep * float(res_s$ID$*res_s$ID$)\n\ + alphaV = kinViscosity_s$ID$ * s$ID$.timestep * float(res_s$ID$*res_s$ID$)\n\ setWallBcs(flags=flags_s$ID$, vel=vel_s$ID$, obvel=None if using_fractions_s$ID$ else obvel_s$ID$, phiObs=phiObs_s$ID$, fractions=fractions_s$ID$)\n\ cgSolveDiffusion(flags_s$ID$, vel_s$ID$, alphaV)\n\ \n\ @@ -315,7 +324,11 @@ def liquid_step_$ID$():\n\ curvature_s$ID$.clamp(-1.0, 1.0)\n\ \n\ setWallBcs(flags=flags_s$ID$, vel=vel_s$ID$, obvel=None if using_fractions_s$ID$ else obvel_s$ID$, phiObs=phiObs_s$ID$, fractions=fractions_s$ID$)\n\ + if using_viscosity_s$ID$:\n\ + viscosity_s$ID$.setConst(viscosityValue_s$ID$)\n\ + applyViscosity(flags=flags_s$ID$, phi=phi_s$ID$, vel=vel_s$ID$, volumes=volumes_s$ID$, viscosity=viscosity_s$ID$)\n\ \n\ + setWallBcs(flags=flags_s$ID$, vel=vel_s$ID$, obvel=None if using_fractions_s$ID$ else obvel_s$ID$, phiObs=phiObs_s$ID$, fractions=fractions_s$ID$)\n\ if using_guiding_s$ID$:\n\ mantaMsg('Guiding and pressure')\n\ PD_fluid_guiding(vel=vel_s$ID$, velT=velT_s$ID$, flags=flags_s$ID$, phi=phi_s$ID$, curv=curvature_s$ID$, surfTens=surfaceTension_s$ID$, fractions=fractions_s$ID$, weight=weightGuide_s$ID$, blurRadius=beta_sg$ID$, pressure=pressure_s$ID$, tau=tau_sg$ID$, sigma=sigma_sg$ID$, theta=theta_sg$ID$, zeroPressureFixing=domainClosed_s$ID$)\n\ diff --git a/release/scripts/startup/bl_ui/properties_physics_fluid.py b/release/scripts/startup/bl_ui/properties_physics_fluid.py index 3afe7a47028..5c3975ffb76 100644 --- a/release/scripts/startup/bl_ui/properties_physics_fluid.py +++ b/release/scripts/startup/bl_ui/properties_physics_fluid.py @@ -1006,6 +1006,46 @@ class PHYSICS_PT_particles(PhysicButtonsPanel, Panel): split.operator("fluid.free_particles", text="Free Particles") +class PHYSICS_PT_viscosity(PhysicButtonsPanel, Panel): + bl_label = "Viscosity" + bl_parent_id = 'PHYSICS_PT_liquid' + bl_options = {'DEFAULT_CLOSED'} + COMPAT_ENGINES = {'BLENDER_RENDER', 'BLENDER_EEVEE', 'BLENDER_WORKBENCH'} + + @classmethod + def poll(cls, context): + # Fluid viscosity only enabled for liquids + if not PhysicButtonsPanel.poll_liquid_domain(context): + return False + + return (context.engine in cls.COMPAT_ENGINES) + + def draw_header(self, context): + md = context.fluid.domain_settings + domain = context.fluid.domain_settings + is_baking_any = domain.is_cache_baking_any + has_baked_any = domain.has_cache_baked_any + self.layout.enabled = not is_baking_any and not has_baked_any + self.layout.prop(md, "use_viscosity", text="") + + def draw(self, context): + layout = self.layout + layout.use_property_split = True + + domain = context.fluid.domain_settings + layout.active = domain.use_viscosity + + is_baking_any = domain.is_cache_baking_any + has_baked_any = domain.has_cache_baked_any + has_baked_data = domain.has_cache_baked_data + + flow = layout.grid_flow(row_major=True, columns=0, even_columns=True, even_rows=False, align=False) + flow.enabled = not is_baking_any and not has_baked_any and not has_baked_data + + col = flow.column(align=True) + col.prop(domain, "viscosity_value", text="Strength") + + class PHYSICS_PT_diffusion(PhysicButtonsPanel, Panel): bl_label = "Diffusion" bl_parent_id = 'PHYSICS_PT_liquid' @@ -1470,6 +1510,7 @@ classes = ( PHYSICS_PT_noise, PHYSICS_PT_fire, PHYSICS_PT_liquid, + PHYSICS_PT_viscosity, PHYSICS_PT_diffusion, PHYSICS_PT_particles, PHYSICS_PT_mesh, diff --git a/source/blender/blenkernel/BKE_blender_version.h b/source/blender/blenkernel/BKE_blender_version.h index 1ed4d1183a1..06a11fc8d1d 100644 --- a/source/blender/blenkernel/BKE_blender_version.h +++ b/source/blender/blenkernel/BKE_blender_version.h @@ -39,7 +39,7 @@ extern "C" { /* Blender file format version. */ #define BLENDER_FILE_VERSION BLENDER_VERSION -#define BLENDER_FILE_SUBVERSION 8 +#define BLENDER_FILE_SUBVERSION 9 /* Minimum Blender version that supports reading file written with the current * version. Older Blender versions will test this and show a warning if the file diff --git a/source/blender/blenkernel/intern/fluid.c b/source/blender/blenkernel/intern/fluid.c index 1bdc2a5dfd2..b93e6148b86 100644 --- a/source/blender/blenkernel/intern/fluid.c +++ b/source/blender/blenkernel/intern/fluid.c @@ -5011,6 +5011,9 @@ void BKE_fluid_modifier_copy(const struct FluidModifierData *fmd, tfds->sys_particle_maximum = fds->sys_particle_maximum; tfds->simulation_method = fds->simulation_method; + /* viscosity options */ + tfds->viscosity_value = fds->viscosity_value; + /* diffusion options*/ tfds->surface_tension = fds->surface_tension; tfds->viscosity_base = fds->viscosity_base; diff --git a/source/blender/blenloader/intern/versioning_290.c b/source/blender/blenloader/intern/versioning_290.c index 8a85278a931..63ef436d8e2 100644 --- a/source/blender/blenloader/intern/versioning_290.c +++ b/source/blender/blenloader/intern/versioning_290.c @@ -1430,18 +1430,7 @@ void blo_do_versions_290(FileData *fd, Library *UNUSED(lib), Main *bmain) } } - /** - * Versioning code until next subversion bump goes here. - * - * \note Be sure to check when bumping the version: - * - "versioning_userdef.c", #blo_do_versions_userdef - * - "versioning_userdef.c", #do_versions_theme - * - * \note Keep this message at the bottom of the function. - */ - { - /* Keep this block, even when empty. */ - + if (!MAIN_VERSION_ATLEAST(bmain, 292, 9)) { FOREACH_NODETREE_BEGIN (bmain, ntree, id) { if (ntree->type == NTREE_GEOMETRY) { LISTBASE_FOREACH (bNode *, node, &ntree->nodes) { @@ -1474,5 +1463,32 @@ void blo_do_versions_290(FileData *fd, Library *UNUSED(lib), Main *bmain) } } } + + /* Ensure that new viscosity strength field is initialized correctly. */ + if (!DNA_struct_elem_find(fd->filesdna, "FluidModifierData", "float", "viscosity_value")) { + LISTBASE_FOREACH (Object *, ob, &bmain->objects) { + LISTBASE_FOREACH (ModifierData *, md, &ob->modifiers) { + if (md->type == eModifierType_Fluid) { + FluidModifierData *fmd = (FluidModifierData *)md; + if (fmd->domain != NULL) { + fmd->domain->viscosity_value = 0.05; + } + } + } + } + } + } + + /** + * Versioning code until next subversion bump goes here. + * + * \note Be sure to check when bumping the version: + * - "versioning_userdef.c", #blo_do_versions_userdef + * - "versioning_userdef.c", #do_versions_theme + * + * \note Keep this message at the bottom of the function. + */ + { + /* Keep this block, even when empty. */ } } diff --git a/source/blender/blenloader/intern/versioning_userdef.c b/source/blender/blenloader/intern/versioning_userdef.c index f1572c563bf..b9f09e0326d 100644 --- a/source/blender/blenloader/intern/versioning_userdef.c +++ b/source/blender/blenloader/intern/versioning_userdef.c @@ -820,6 +820,12 @@ void blo_do_versions_userdef(UserDef *userdef) userdef->uiflag &= ~USER_UIFLAG_UNUSED_3; } + if (!USER_VERSION_ATLEAST(292, 9)) { + if (BLI_listbase_is_empty(&userdef->asset_libraries)) { + BKE_preferences_asset_library_default_add(userdef); + } + } + /** * Versioning code until next subversion bump goes here. * @@ -831,9 +837,6 @@ void blo_do_versions_userdef(UserDef *userdef) */ { /* Keep this block, even when empty. */ - if (BLI_listbase_is_empty(&userdef->asset_libraries)) { - BKE_preferences_asset_library_default_add(userdef); - } } LISTBASE_FOREACH (bTheme *, btheme, &userdef->themes) { diff --git a/source/blender/makesdna/DNA_fluid_defaults.h b/source/blender/makesdna/DNA_fluid_defaults.h index 2ee83cf3387..9454342654c 100644 --- a/source/blender/makesdna/DNA_fluid_defaults.h +++ b/source/blender/makesdna/DNA_fluid_defaults.h @@ -113,6 +113,7 @@ .flip_ratio = 0.97f, \ .sys_particle_maximum = 0, \ .simulation_method = FLUID_DOMAIN_METHOD_FLIP, \ + .viscosity_value = 0.05f, \ .surface_tension = 0.0f, \ .viscosity_base = 1.0f, \ .viscosity_exponent = 6.0f, \ diff --git a/source/blender/makesdna/DNA_fluid_types.h b/source/blender/makesdna/DNA_fluid_types.h index 98007a2b0b1..2ba0d15fbb3 100644 --- a/source/blender/makesdna/DNA_fluid_types.h +++ b/source/blender/makesdna/DNA_fluid_types.h @@ -52,6 +52,7 @@ enum { FLUID_DOMAIN_DELETE_IN_OBSTACLE = (1 << 14), /* Delete fluid inside obstacles. */ FLUID_DOMAIN_USE_DIFFUSION = (1 << 15), /* Use diffusion (e.g. viscosity, surface tension). */ FLUID_DOMAIN_USE_RESUMABLE_CACHE = (1 << 16), /* Determine if cache should be resumable. */ + FLUID_DOMAIN_USE_VISCOSITY = (1 << 17), /* Use viscosity. */ }; /** @@ -289,7 +290,7 @@ enum { #define FLUID_NAME_GUIDING "fluid_guiding" /* Fluid object names.*/ -#define FLUID_NAME_FLAGS "flags" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_FLAGS "flags" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_VELOCITY "velocity" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_VEL "vel" #define FLUID_NAME_VELOCITYTMP "velocity_previous" /* == OpenVDB grid attribute name. */ @@ -300,7 +301,7 @@ enum { #define FLUID_NAME_PHIOBS "phi_obstacle" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_PHISIN "phiSIn" #define FLUID_NAME_PHIIN "phi_inflow" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_PHIOUT "phi_out" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_PHIOUT "phi_out" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_FORCES "forces" #define FLUID_NAME_FORCE_X "x_force" #define FLUID_NAME_FORCE_Y "y_force" @@ -322,37 +323,37 @@ enum { #define FLUID_NAME_PHIOUTIN "phi_out_inflow" /* Smoke object names. */ -#define FLUID_NAME_SHADOW "shadow" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_SHADOW "shadow" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_EMISSION "emission" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_EMISSIONIN "emissionIn" -#define FLUID_NAME_DENSITY "density" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_DENSITY "density" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_DENSITYIN "density_inflow" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_HEAT "heat" #define FLUID_NAME_HEATIN "heatIn" -#define FLUID_NAME_TEMPERATURE "temperature" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_TEMPERATURE "temperature" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_TEMPERATUREIN "temperature_inflow" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_COLORR "color_r" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_COLORG "color_g" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_COLORB "color_b" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_COLORRIN "color_r_inflow" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_COLORGIN "color_g_inflow" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_COLORBIN "color_b_inflow" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_FLAME "flame" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_FUEL "fuel" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_REACT "react" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_FUELIN "fuel_inflow" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_REACTIN "react_inflow" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_COLORR "color_r" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_COLORG "color_g" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_COLORB "color_b" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_COLORRIN "color_r_inflow" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_COLORGIN "color_g_inflow" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_COLORBIN "color_b_inflow" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_FLAME "flame" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_FUEL "fuel" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_REACT "react" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_FUELIN "fuel_inflow" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_REACTIN "react_inflow" /* == OpenVDB grid attribute name. */ /* Liquid object names. */ #define FLUID_NAME_PHIPARTS "phi_particles" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_PHI "phi" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_PHITMP "phi_previous" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_PHI "phi" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_PHITMP "phi_previous" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_VELOCITYOLD "velOld" #define FLUID_NAME_VELOCITYPARTS "velParts" #define FLUID_NAME_MAPWEIGHTS "mapWeights" #define FLUID_NAME_PP "pp" #define FLUID_NAME_PVEL "pVel" -#define FLUID_NAME_PARTS "particles" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_PARTS "particles" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_PARTSVELOCITY "particles_velocity" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_PINDEX "pindex" #define FLUID_NAME_GPI "gpi" @@ -375,8 +376,8 @@ enum { #define FLUID_NAME_TEXTURE_U2 "textureU2" #define FLUID_NAME_TEXTURE_V2 "textureV2" #define FLUID_NAME_TEXTURE_W2 "textureW2" -#define FLUID_NAME_UV0 "uv_grid_0" /* == OpenVDB grid attribute name. */ -#define FLUID_NAME_UV1 "uv_grid_1" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_UV0 "uv_grid_0" /* == OpenVDB grid attribute name. */ +#define FLUID_NAME_UV1 "uv_grid_1" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_COLORR_NOISE "color_r_noise" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_COLORG_NOISE "color_g_noise" /* == OpenVDB grid attribute name. */ #define FLUID_NAME_COLORB_NOISE "color_b_noise" /* == OpenVDB grid attribute name. */ @@ -596,6 +597,10 @@ typedef struct FluidDomainSettings { short simulation_method; char _pad4[6]; + /* Viscosity options. */ + float viscosity_value; + char _pad5[4]; + /* Diffusion options. */ float surface_tension; float viscosity_base; @@ -610,7 +615,7 @@ typedef struct FluidDomainSettings { int mesh_scale; int totvert; short mesh_generator; - char _pad5[6]; /* Unused. */ + char _pad6[6]; /* Unused. */ /* Secondary particle options. */ int particle_type; @@ -631,7 +636,7 @@ typedef struct FluidDomainSettings { int sndparticle_update_radius; char sndparticle_boundary; char sndparticle_combined_export; - char _pad6[6]; /* Unused. */ + char _pad7[6]; /* Unused. */ /* Fluid guiding options. */ float guide_alpha; /* Guiding weight scalar (determines strength). */ @@ -639,7 +644,7 @@ typedef struct FluidDomainSettings { float guide_vel_factor; /* Multiply guiding velocity by this factor. */ int guide_res[3]; /* Res for velocity guide grids - independent from base res. */ short guide_source; - char _pad7[2]; /* Unused. */ + char _pad8[2]; /* Unused. */ /* Cache options. */ int cache_frame_start; @@ -659,7 +664,7 @@ typedef struct FluidDomainSettings { char error[64]; /* Bake error description. */ short cache_type; char cache_id[4]; /* Run-time only */ - char _pad8[2]; + char _pad9[2]; /* Unused. */ /* Time options. */ float dt; @@ -694,19 +699,19 @@ typedef struct FluidDomainSettings { char interp_method; char gridlines_color_field; /* Simulation field used to color map onto gridlines. */ char gridlines_cell_filter; - char _pad9[7]; + char _pad10[7]; /* Unused. */ /* OpenVDB cache options. */ int openvdb_compression; float clipping; char openvdb_data_depth; - char _pad10[7]; /* Unused. */ + char _pad11[7]; /* Unused. */ /* -- Deprecated / unsed options (below). -- */ /* View options. */ int viewsettings; - char _pad11[4]; /* Unused. */ + char _pad12[4]; /* Unused. */ /* Pointcache options. */ /* Smoke uses only one cache from now on (index [0]), but keeping the array for now for reading @@ -716,7 +721,7 @@ typedef struct FluidDomainSettings { int cache_comp; int cache_high_comp; char cache_file_format; - char _pad12[7]; /* Unused. */ + char _pad13[7]; /* Unused. */ } FluidDomainSettings; diff --git a/source/blender/makesrna/intern/rna_fluid.c b/source/blender/makesrna/intern/rna_fluid.c index bb8280ede91..afe564dff0a 100644 --- a/source/blender/makesrna/intern/rna_fluid.c +++ b/source/blender/makesrna/intern/rna_fluid.c @@ -1929,6 +1929,22 @@ static void rna_def_fluid_domain_settings(BlenderRNA *brna) "Maximum number of fluid particles that are allowed in this simulation"); RNA_def_property_update(prop, NC_OBJECT | ND_MODIFIER, "rna_Fluid_datacache_reset"); + /* viscosity options */ + + prop = RNA_def_property(srna, "use_viscosity", PROP_BOOLEAN, PROP_NONE); + RNA_def_property_boolean_sdna(prop, NULL, "flags", FLUID_DOMAIN_USE_VISCOSITY); + RNA_def_property_ui_text(prop, "Use Viscosity", "Enable fluid viscosity settings"); + RNA_def_property_update(prop, NC_OBJECT | ND_MODIFIER, "rna_Fluid_datacache_reset"); + + prop = RNA_def_property(srna, "viscosity_value", PROP_FLOAT, PROP_NONE); + RNA_def_property_range(prop, 0.0, 10.0); + RNA_def_property_ui_range(prop, 0.0, 5.0, 0.01, 3); + RNA_def_property_ui_text(prop, + "Strength", + "Viscosity of liquid (higher values result in more viscous fluids, a " + "value of 0 will still apply some viscosity)"); + RNA_def_property_update(prop, NC_OBJECT | ND_MODIFIER, "rna_Fluid_datacache_reset"); + /* diffusion options */ prop = RNA_def_property(srna, "use_diffusion", PROP_BOOLEAN, PROP_NONE); |