Fluid: Added new viscosity solver

Mainly updated the Mantaflow version. It includes the new viscosity solver plugin based on the method from 'Accurate Viscous Free Surfaces for Buckling, Coiling, and Rotating Liquids' (Batty & Bridson).

In the UI, this update adds a new 'Viscosity' section to the fluid modifier UI (liquid domains only). For now, there is a single 'strength' value to control the viscosity of liquids.

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

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

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,
+                         std::vector<Grid<Real> *> matrixAVec,
-                         Grid<Real> *pAi,
+                         std::vector<Grid<Real> *> rhsVec)
-                         Grid<Real> *pAj,
-                         Grid<Real> *pAk)
     : GridCgInterface(),
       mInited(false),
       mIterations(0),
@@ -419,10 +417,8 @@ GridCg<APPLYMAT>::GridCg(Grid<Real> &dst,
       mSearch(search),
       mFlags(flags),
       mTmp(tmp),
-      mpA0(pA0),
+      mMatrixA(matrixAVec),
-      mpAi(pAi),
+      mVecRhs(rhsVec),
-      mpAj(pAj),
-      mpAk(pAk),
       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(
+    InitPreconditionIncompCholesky(mFlags,
-        mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
+                                   *mpPCA0,
-    ApplyPreconditionIncompCholesky(
+                                   *mpPCAi,
-        mTmp, mResidual, mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
+                                   *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(
+    ApplyPreconditionIncompCholesky(mTmp,
-        mTmp, mResidual, mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
+                                    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())
+    vector<Grid<Real> *> matA{&A0, &Ai, &Aj};
-      gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+
-    else
+    if (flags.is3D()) {
-      gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+      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())
+    if (flags.is3D()) {
-      gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+      matA.push_back(&Ak);
-    else
+      gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, matA);
-      gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+    }
+    else {
+      gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, matA);
+    }
+
     gcg->setAccuracy(cgAccuracy);
 
     // diffuse every component separately

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,
+         std::vector<Grid<Real> *> matrixAVec,
-         Grid<Real> *pAi,
+         std::vector<Grid<Real> *> rhsVec = {});
-         Grid<Real> *pAj,
+  ~GridCg(){};
-         Grid<Real> *pAk);
-  ~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,
+              const std::vector<Grid<Real> *> matrixA,
-              Grid<Real> &Ai,
+              const std::vector<Grid<Real> *> vecRhs)
-              Grid<Real> &Aj,
+      : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
-              Grid<Real> &Ak)
-      : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), A0(A0), Ai(Ai), Aj(Aj), Ak(Ak)
   {
     runMessage();
     run();
@@ -167,11 +164,18 @@ struct ApplyMatrix : public KernelBase {
                  const FlagGrid &flags,
                  Grid<Real> &dst,
                  const Grid<Real> &src,
-                 Grid<Real> &A0,
+                 const std::vector<Grid<Real> *> matrixA,
-                 Grid<Real> &Ai,
+                 const std::vector<Grid<Real> *> vecRhs) const
-                 Grid<Real> &Aj,
-                 Grid<Real> &Ak) 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()
+  inline const std::vector<Grid<Real> *> &getArg3()
   {
-    return A0;
+    return matrixA;
   }
-  typedef Grid<Real> type3;
+  typedef std::vector<Grid<Real> *> type3;
-  inline Grid<Real> &getArg4()
+  inline const std::vector<Grid<Real> *> &getArg4()
   {
-    return Ai;
+    return vecRhs;
   }
-  typedef Grid<Real> type4;
+  typedef std::vector<Grid<Real> *> type4;
-  inline Grid<Real> &getArg5()
-  {
-    return Aj;
-  }
-  typedef Grid<Real> type5;
-  inline Grid<Real> &getArg6()
-  {
-    return Ak;
-  }
-  typedef Grid<Real> type6;
   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;
+  const std::vector<Grid<Real> *> matrixA;
-  Grid<Real> &Ai;
+  const std::vector<Grid<Real> *> vecRhs;
-  Grid<Real> &Aj;
-  Grid<Real> &Ak;
 };
 
 //! 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,
+                const std::vector<Grid<Real> *> matrixA,
-                Grid<Real> &Ai,
+                const std::vector<Grid<Real> *> vecRhs)
-                Grid<Real> &Aj,
+      : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
-                Grid<Real> &Ak)
-      : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), A0(A0), Ai(Ai), Aj(Aj), Ak(Ak)
   {
     runMessage();
     run();
@@ -260,12 +250,16 @@ struct ApplyMatrix2D : public KernelBase {
                  const FlagGrid &flags,
                  Grid<Real> &dst,
                  const Grid<Real> &src,
-                 Grid<Real> &A0,
+                 const std::vector<Grid<Real> *> matrixA,
-                 Grid<Real> &Ai,
+                 const std::vector<Grid<Real> *> vecRhs) const
-                 Grid<Real> &Aj,
-                 Grid<Real> &Ak) 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,26 +284,16 @@ 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;
+  typedef std::vector<Grid<Real> *> type3;
-  inline Grid<Real> &getArg4()
+  inline const std::vector<Grid<Real> *> &getArg4()
   {
-    return Ai;
+    return vecRhs;
   }
-  typedef Grid<Real> type4;
+  typedef std::vector<Grid<Real> *> type4;
-  inline Grid<Real> &getArg5()
-  {
-    return Aj;
-  }
-  typedef Grid<Real> type5;
-  inline Grid<Real> &getArg6()
-  {
-    return Ak;
-  }
-  typedef Grid<Real> type6;
   void runMessage()
   {
     debMsg("Executing kernel ApplyMatrix2D ", 3);
@@ -320,7 +304,7 @@ struct ApplyMatrix2D : 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()
   {
@@ -329,12 +313,358 @@ struct ApplyMatrix2D : public KernelBase {
   const FlagGrid &flags;
   Grid<Real> &dst;
   const Grid<Real> &src;
-  Grid<Real> &A0;
+  const std::vector<Grid<Real> *> matrixA;
-  Grid<Real> &Ai;
+  const std::vector<Grid<Real> *> vecRhs;
-  Grid<Real> &Aj;
-  Grid<Real> &Ak;
 };
 
+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();
+  }
+  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 ApplyMatrixViscosityU ", 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 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 {

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

extern/mantaflow/preprocessed/gitinfo.h

@@ -1,3 +1,3 @@
 
 
-#define MANTA_GIT_VERSION "commit 327917cd59b03bef3a953b5f58fc1637b3a83e01"
+#define MANTA_GIT_VERSION "commit e2285cb9bc492987f728123be6cfc1fe11fe73d6"

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);
+      minVal = min(vel(i, j, k).x, forceMACX);
-      Real max = std::max(vel(i, j, k).x, forceMACX);
+      maxVal = max(vel(i, j, k).x, forceMACX);
-      Real sum = vel(i, j, k).x + forceMACX;
+      sum = vel(i, j, k).x + forceMACX;
-      vel(i, j, k).x = (forceMACX > 0) ? std::min(sum, max) : std::max(sum, min);
+      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);
+      minVal = min(vel(i, j, k).y, forceMACY);
-      Real max = std::max(vel(i, j, k).y, forceMACY);
+      maxVal = max(vel(i, j, k).y, forceMACY);
-      Real sum = vel(i, j, k).y + forceMACY;
+      sum = vel(i, j, k).y + forceMACY;
-      vel(i, j, k).y = (forceMACY > 0) ? std::min(sum, max) : std::max(sum, min);
+      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);
+      minVal = min(vel(i, j, k).z, forceMACZ);
-      Real max = std::max(vel(i, j, k).z, forceMACZ);
+      maxVal = max(vel(i, j, k).z, forceMACZ);
-      Real sum = vel(i, j, k).z + forceMACZ;
+      sum = vel(i, j, k).z + forceMACZ;
-      vel(i, j, k).z = (forceMACZ > 0) ? std::min(sum, max) : std::max(sum, min);
+      vel(i, j, k).z = (forceMACZ > 0) ? min(sum, maxVal) : max(sum, minVal);
     }
   }
   inline const FlagGrid &getArg0()

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())
+  vector<Grid<Real> *> matA{&A0, &Ai, &Aj};
-    gcg = new GridCg<ApplyMatrix>(pressure, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+
-  else
+  if (vel.is3D()) {
-    gcg = new GridCg<ApplyMatrix2D>(
+    matA.push_back(&Ak);
-        pressure, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &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);

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(

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())
+  vector<Grid<Real> *> matA{&A0, &Ai, &Aj};
-    gcg = new GridCg<ApplyMatrix>(out, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+
-  else
+  if (flags.is3D()) {
-    gcg = new GridCg<ApplyMatrix2D>(out, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+    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);
 

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();

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;

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;

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";

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\

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,

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

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;

source/blender/blenloader/intern/versioning_290.c

@@ -1430,18 +1430,7 @@ void blo_do_versions_290(FileData *fd, Library *UNUSED(lib), Main *bmain)
     }
   }
 
-  /**
+  if (!MAIN_VERSION_ATLEAST(bmain, 292, 9)) {
-   * 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. */
-
     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. */
   }
 }

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) {

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, \

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_COLORR "color_r" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_COLORG "color_g"                   /* == 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_COLORB "color_b" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_COLORRIN "color_r_inflow"          /* == 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_COLORGIN "color_g_inflow" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_COLORBIN "color_b_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_FLAME "flame" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_FUEL "fuel"                        /* == 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_REACT "react" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_FUELIN "fuel_inflow"               /* == 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_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_PHI "phi" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_PHITMP "phi_previous"    /* == 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_UV0 "uv_grid_0" /* == OpenVDB grid attribute name. */
-#define FLUID_NAME_UV1 "uv_grid_1"              /* == 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;
 

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);