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.

12 files changed, 1965 insertions, 123 deletions

diff --git a/extern/mantaflow/CMakeLists.txt b/extern/mantaflow/CMakeLists.txt
index ccf272650e3..82bf95a9742 100644
--- a/extern/mantaflow/CMakeLists.txt
+++ b/extern/mantaflow/CMakeLists.txt
@@ -200,6 +200,7 @@ set(SRC
   ${MANTA_PP}/plugin/ptsplugins.cpp
   ${MANTA_PP}/plugin/secondaryparticles.cpp
   ${MANTA_PP}/plugin/surfaceturbulence.cpp
+  ${MANTA_PP}/plugin/viscosity.cpp
   ${MANTA_PP}/plugin/vortexplugins.cpp
   ${MANTA_PP}/plugin/waveletturbulence.cpp
   ${MANTA_PP}/plugin/waves.cpp
diff --git a/extern/mantaflow/helper/util/rcmatrix.h b/extern/mantaflow/helper/util/rcmatrix.h
index f1f0efe6416..330fd1f64f7 100644
--- a/extern/mantaflow/helper/util/rcmatrix.h
+++ b/extern/mantaflow/helper/util/rcmatrix.h
@@ -1035,7 +1035,7 @@ template<class N, class T> struct RCFixedMatrix {
 typedef RCMatrix<int, Real> Matrix;
 typedef RCFixedMatrix<int, Real> FixedMatrix;
 
-}  // namespace Manta
+}
 
 #undef parallel_for
 #undef parallel_end
diff --git a/extern/mantaflow/preprocessed/conjugategrad.cpp b/extern/mantaflow/preprocessed/conjugategrad.cpp
index da82a412a97..bdcceb29520 100644
--- a/extern/mantaflow/preprocessed/conjugategrad.cpp
+++ b/extern/mantaflow/preprocessed/conjugategrad.cpp
@@ -397,7 +397,7 @@ struct UpdateSearchVec : public KernelBase {
 };
 
 //*****************************************************************************
-//  CG class
+// CG class
 
 template<class APPLYMAT>
 GridCg<APPLYMAT>::GridCg(Grid<Real> &dst,
@@ -406,10 +406,8 @@ GridCg<APPLYMAT>::GridCg(Grid<Real> &dst,
                          Grid<Real> &search,
                          const FlagGrid &flags,
                          Grid<Real> &tmp,
-                         Grid<Real> *pA0,
-                         Grid<Real> *pAi,
-                         Grid<Real> *pAj,
-                         Grid<Real> *pAk)
+                         std::vector<Grid<Real> *> matrixAVec,
+                         std::vector<Grid<Real> *> rhsVec)
     : GridCgInterface(),
       mInited(false),
       mIterations(0),
@@ -419,10 +417,8 @@ GridCg<APPLYMAT>::GridCg(Grid<Real> &dst,
       mSearch(search),
       mFlags(flags),
       mTmp(tmp),
-      mpA0(pA0),
-      mpAi(pAi),
-      mpAj(pAj),
-      mpAk(pAk),
+      mMatrixA(matrixAVec),
+      mVecRhs(rhsVec),
       mPcMethod(PC_None),
       mpPCA0(nullptr),
       mpPCAi(nullptr),
@@ -445,19 +441,37 @@ template<class APPLYMAT> void GridCg<APPLYMAT>::doInit()
 
   if (mPcMethod == PC_ICP) {
     assertMsg(mDst.is3D(), "ICP only supports 3D grids so far");
-    InitPreconditionIncompCholesky(
-        mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
-    ApplyPreconditionIncompCholesky(
-        mTmp, mResidual, mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
+    InitPreconditionIncompCholesky(mFlags,
+                                   *mpPCA0,
+                                   *mpPCAi,
+                                   *mpPCAj,
+                                   *mpPCAk,
+                                   *mMatrixA[0],
+                                   *mMatrixA[1],
+                                   *mMatrixA[2],
+                                   *mMatrixA[3]);
+    ApplyPreconditionIncompCholesky(mTmp,
+                                    mResidual,
+                                    mFlags,
+                                    *mpPCA0,
+                                    *mpPCAi,
+                                    *mpPCAj,
+                                    *mpPCAk,
+                                    *mMatrixA[0],
+                                    *mMatrixA[1],
+                                    *mMatrixA[2],
+                                    *mMatrixA[3]);
   }
   else if (mPcMethod == PC_mICP) {
     assertMsg(mDst.is3D(), "mICP only supports 3D grids so far");
-    InitPreconditionModifiedIncompCholesky2(mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk);
+    InitPreconditionModifiedIncompCholesky2(
+        mFlags, *mpPCA0, *mMatrixA[0], *mMatrixA[1], *mMatrixA[2], *mMatrixA[3]);
     ApplyPreconditionModifiedIncompCholesky2(
-        mTmp, mResidual, mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk);
+        mTmp, mResidual, mFlags, *mpPCA0, *mMatrixA[0], *mMatrixA[1], *mMatrixA[2], *mMatrixA[3]);
   }
   else if (mPcMethod == PC_MGP) {
-    InitPreconditionMultigrid(mMG, *mpA0, *mpAi, *mpAj, *mpAk, mAccuracy);
+    InitPreconditionMultigrid(
+        mMG, *mMatrixA[0], *mMatrixA[1], *mMatrixA[2], *mMatrixA[3], mAccuracy);
     ApplyPreconditionMultigrid(mMG, mTmp, mResidual);
   }
   else {
@@ -465,7 +479,6 @@ template<class APPLYMAT> void GridCg<APPLYMAT>::doInit()
   }
 
   mSearch.copyFrom(mTmp);
-
   mSigma = GridDotProduct(mTmp, mResidual);
 }
 
@@ -480,7 +493,7 @@ template<class APPLYMAT> bool GridCg<APPLYMAT>::iterate()
   // this could reinterpret the mpA pointers (not so clean right now)
   // tmp = applyMat(search)
 
-  APPLYMAT(mFlags, mTmp, mSearch, *mpA0, *mpAi, *mpAj, *mpAk);
+  APPLYMAT(mFlags, mTmp, mSearch, mMatrixA, mVecRhs);
 
   // alpha = sigma/dot(tmp, search)
   Real dp = GridDotProduct(mTmp, mSearch);
@@ -492,11 +505,20 @@ template<class APPLYMAT> bool GridCg<APPLYMAT>::iterate()
   gridScaledAdd<Real, Real>(mResidual, mTmp, -alpha);  // residual += tmp * -alpha
 
   if (mPcMethod == PC_ICP)
-    ApplyPreconditionIncompCholesky(
-        mTmp, mResidual, mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
+    ApplyPreconditionIncompCholesky(mTmp,
+                                    mResidual,
+                                    mFlags,
+                                    *mpPCA0,
+                                    *mpPCAi,
+                                    *mpPCAj,
+                                    *mpPCAk,
+                                    *mMatrixA[0],
+                                    *mMatrixA[1],
+                                    *mMatrixA[2],
+                                    *mMatrixA[3]);
   else if (mPcMethod == PC_mICP)
     ApplyPreconditionModifiedIncompCholesky2(
-        mTmp, mResidual, mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk);
+        mTmp, mResidual, mFlags, *mpPCA0, *mMatrixA[0], *mMatrixA[1], *mMatrixA[2], *mMatrixA[3]);
   else if (mPcMethod == PC_MGP)
     ApplyPreconditionMultigrid(mMG, mTmp, mResidual);
   else
@@ -584,13 +606,15 @@ void GridCg<APPLYMAT>::setMGPreconditioner(PreconditionType method, GridMg *MG)
   assertMsg(method == PC_MGP, "GridCg<APPLYMAT>::setMGPreconditioner: Invalid method specified.");
 
   mPcMethod = method;
-
   mMG = MG;
 }
 
 // explicit instantiation
 template class GridCg<ApplyMatrix>;
 template class GridCg<ApplyMatrix2D>;
+template class GridCg<ApplyMatrixViscosityU>;
+template class GridCg<ApplyMatrixViscosityV>;
+template class GridCg<ApplyMatrixViscosityW>;
 
 //*****************************************************************************
 // diffusion for real and vec grids, e.g. for viscosity
@@ -638,10 +662,15 @@ void cgSolveDiffusion(const FlagGrid &flags,
   if (grid.getType() & GridBase::TypeReal) {
     Grid<Real> &u = ((Grid<Real> &)grid);
     rhs.copyFrom(u);
-    if (flags.is3D())
-      gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
-    else
-      gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+    vector<Grid<Real> *> matA{&A0, &Ai, &Aj};
+
+    if (flags.is3D()) {
+      matA.push_back(&Ak);
+      gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, matA);
+    }
+    else {
+      gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, matA);
+    }
 
     gcg->setAccuracy(cgAccuracy);
     gcg->solve(maxIter);
@@ -653,12 +682,17 @@ void cgSolveDiffusion(const FlagGrid &flags,
   else if ((grid.getType() & GridBase::TypeVec3) || (grid.getType() & GridBase::TypeMAC)) {
     Grid<Vec3> &vec = ((Grid<Vec3> &)grid);
     Grid<Real> u(parent);
+    vector<Grid<Real> *> matA{&A0, &Ai, &Aj};
 
     // core solve is same as for a regular real grid
-    if (flags.is3D())
-      gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
-    else
-      gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+    if (flags.is3D()) {
+      matA.push_back(&Ak);
+      gcg = new GridCg<ApplyMatrix>(u, rhs, residual, search, flags, tmp, matA);
+    }
+    else {
+      gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, matA);
+    }
+
     gcg->setAccuracy(cgAccuracy);
 
     // diffuse every component separately
diff --git a/extern/mantaflow/preprocessed/conjugategrad.h b/extern/mantaflow/preprocessed/conjugategrad.h
index 58ccff28179..4974aa6a4d6 100644
--- a/extern/mantaflow/preprocessed/conjugategrad.h
+++ b/extern/mantaflow/preprocessed/conjugategrad.h
@@ -78,13 +78,9 @@ template<class APPLYMAT> class GridCg : public GridCgInterface {
          Grid<Real> &search,
          const FlagGrid &flags,
          Grid<Real> &tmp,
-         Grid<Real> *A0,
-         Grid<Real> *pAi,
-         Grid<Real> *pAj,
-         Grid<Real> *pAk);
-  ~GridCg()
-  {
-  }
+         std::vector<Grid<Real> *> matrixAVec,
+         std::vector<Grid<Real> *> rhsVec = {});
+  ~GridCg(){};
 
   void doInit();
   bool iterate();
@@ -133,7 +129,10 @@ template<class APPLYMAT> class GridCg : public GridCgInterface {
   const FlagGrid &mFlags;
   Grid<Real> &mTmp;
 
-  Grid<Real> *mpA0, *mpAi, *mpAj, *mpAk;
+  //! shape of A matrix defined here (e.g. diagonal,  positive neighbor cells, etc)
+  std::vector<Grid<Real> *> mMatrixA;
+  //! shape of rhs vector defined here (e.g. 1 rhs for regular fluids solve, 3 rhs for viscosity)
+  std::vector<Grid<Real> *> mVecRhs;
 
   PreconditionType mPcMethod;
   //! preconditioning grids
@@ -154,11 +153,9 @@ struct ApplyMatrix : public KernelBase {
   ApplyMatrix(const FlagGrid &flags,
               Grid<Real> &dst,
               const Grid<Real> &src,
-              Grid<Real> &A0,
-              Grid<Real> &Ai,
-              Grid<Real> &Aj,
-              Grid<Real> &Ak)
-      : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), A0(A0), Ai(Ai), Aj(Aj), Ak(Ak)
+              const std::vector<Grid<Real> *> matrixA,
+              const std::vector<Grid<Real> *> vecRhs)
+      : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
   {
     runMessage();
     run();
@@ -167,11 +164,18 @@ struct ApplyMatrix : public KernelBase {
                  const FlagGrid &flags,
                  Grid<Real> &dst,
                  const Grid<Real> &src,
-                 Grid<Real> &A0,
-                 Grid<Real> &Ai,
-                 Grid<Real> &Aj,
-                 Grid<Real> &Ak) const
+                 const std::vector<Grid<Real> *> matrixA,
+                 const std::vector<Grid<Real> *> vecRhs) const
   {
+    unusedParameter(vecRhs);  // Not needed in this matrix application
+
+    if (matrixA.size() != 4)
+      errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step");
+    Grid<Real> &A0 = *matrixA[0];
+    Grid<Real> &Ai = *matrixA[1];
+    Grid<Real> &Aj = *matrixA[2];
+    Grid<Real> &Ak = *matrixA[3];
+
     if (!flags.isFluid(idx)) {
       dst[idx] = src[idx];
       return;
@@ -196,26 +200,16 @@ struct ApplyMatrix : public KernelBase {
     return src;
   }
   typedef Grid<Real> type2;
-  inline Grid<Real> &getArg3()
-  {
-    return A0;
-  }
-  typedef Grid<Real> type3;
-  inline Grid<Real> &getArg4()
-  {
-    return Ai;
-  }
-  typedef Grid<Real> type4;
-  inline Grid<Real> &getArg5()
+  inline const std::vector<Grid<Real> *> &getArg3()
   {
-    return Aj;
+    return matrixA;
   }
-  typedef Grid<Real> type5;
-  inline Grid<Real> &getArg6()
+  typedef std::vector<Grid<Real> *> type3;
+  inline const std::vector<Grid<Real> *> &getArg4()
   {
-    return Ak;
+    return vecRhs;
   }
-  typedef Grid<Real> type6;
+  typedef std::vector<Grid<Real> *> type4;
   void runMessage()
   {
     debMsg("Executing kernel ApplyMatrix ", 3);
@@ -226,7 +220,7 @@ struct ApplyMatrix : public KernelBase {
   void operator()(const tbb::blocked_range<IndexInt> &__r) const
   {
     for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
-      op(idx, flags, dst, src, A0, Ai, Aj, Ak);
+      op(idx, flags, dst, src, matrixA, vecRhs);
   }
   void run()
   {
@@ -235,10 +229,8 @@ struct ApplyMatrix : public KernelBase {
   const FlagGrid &flags;
   Grid<Real> &dst;
   const Grid<Real> &src;
-  Grid<Real> &A0;
-  Grid<Real> &Ai;
-  Grid<Real> &Aj;
-  Grid<Real> &Ak;
+  const std::vector<Grid<Real> *> matrixA;
+  const std::vector<Grid<Real> *> vecRhs;
 };
 
 //! Kernel: Apply symmetric stored Matrix. 2D version
@@ -247,11 +239,9 @@ struct ApplyMatrix2D : public KernelBase {
   ApplyMatrix2D(const FlagGrid &flags,
                 Grid<Real> &dst,
                 const Grid<Real> &src,
-                Grid<Real> &A0,
-                Grid<Real> &Ai,
-                Grid<Real> &Aj,
-                Grid<Real> &Ak)
-      : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), A0(A0), Ai(Ai), Aj(Aj), Ak(Ak)
+                const std::vector<Grid<Real> *> matrixA,
+                const std::vector<Grid<Real> *> vecRhs)
+      : KernelBase(&flags, 0), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
   {
     runMessage();
     run();
@@ -260,12 +250,16 @@ struct ApplyMatrix2D : public KernelBase {
                  const FlagGrid &flags,
                  Grid<Real> &dst,
                  const Grid<Real> &src,
-                 Grid<Real> &A0,
-                 Grid<Real> &Ai,
-                 Grid<Real> &Aj,
-                 Grid<Real> &Ak) const
+                 const std::vector<Grid<Real> *> matrixA,
+                 const std::vector<Grid<Real> *> vecRhs) const
   {
-    unusedParameter(Ak);  // only there for parameter compatibility with ApplyMatrix
+    unusedParameter(vecRhs);  // Not needed in this matrix application
+
+    if (matrixA.size() != 3)
+      errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step");
+    Grid<Real> &A0 = *matrixA[0];
+    Grid<Real> &Ai = *matrixA[1];
+    Grid<Real> &Aj = *matrixA[2];
 
     if (!flags.isFluid(idx)) {
       dst[idx] = src[idx];
@@ -290,51 +284,387 @@ struct ApplyMatrix2D : public KernelBase {
     return src;
   }
   typedef Grid<Real> type2;
-  inline Grid<Real> &getArg3()
+  inline const std::vector<Grid<Real> *> &getArg3()
   {
-    return A0;
+    return matrixA;
   }
-  typedef Grid<Real> type3;
-  inline Grid<Real> &getArg4()
+  typedef std::vector<Grid<Real> *> type3;
+  inline const std::vector<Grid<Real> *> &getArg4()
   {
-    return Ai;
+    return vecRhs;
   }
-  typedef Grid<Real> type4;
-  inline Grid<Real> &getArg5()
+  typedef std::vector<Grid<Real> *> type4;
+  void runMessage()
   {
-    return Aj;
+    debMsg("Executing kernel ApplyMatrix2D ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r) const
+  {
+    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
+      op(idx, flags, dst, src, matrixA, vecRhs);
   }
-  typedef Grid<Real> type5;
-  inline Grid<Real> &getArg6()
+  void run()
   {
-    return Ak;
+    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
+  }
+  const FlagGrid &flags;
+  Grid<Real> &dst;
+  const Grid<Real> &src;
+  const std::vector<Grid<Real> *> matrixA;
+  const std::vector<Grid<Real> *> vecRhs;
+};
+
+struct ApplyMatrixViscosityU : public KernelBase {
+  ApplyMatrixViscosityU(const FlagGrid &flags,
+                        Grid<Real> &dst,
+                        const Grid<Real> &src,
+                        const std::vector<Grid<Real> *> matrixA,
+                        const std::vector<Grid<Real> *> vecRhs)
+      : KernelBase(&flags, 1), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
+  {
+    runMessage();
+    run();
   }
-  typedef Grid<Real> type6;
+  inline void op(int i,
+                 int j,
+                 int k,
+                 const FlagGrid &flags,
+                 Grid<Real> &dst,
+                 const Grid<Real> &src,
+                 const std::vector<Grid<Real> *> matrixA,
+                 const std::vector<Grid<Real> *> vecRhs) const
+  {
+    if (matrixA.size() != 15)
+      errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step");
+    Grid<Real> &A0 = *matrixA[0];
+    Grid<Real> &Aplusi = *matrixA[1];
+    Grid<Real> &Aplusj = *matrixA[2];
+    Grid<Real> &Aplusk = *matrixA[3];
+    Grid<Real> &Aminusi = *matrixA[4];
+    Grid<Real> &Aminusj = *matrixA[5];
+    Grid<Real> &Aminusk = *matrixA[6];
+
+    if (vecRhs.size() != 2)
+      errMsg("ConjugatedGrad: Invalid rhs vector in apply matrix step");
+    Grid<Real> &srcV = *vecRhs[0];
+    Grid<Real> &srcW = *vecRhs[1];
+
+    dst(i, j, k) = src(i, j, k) * A0(i, j, k) + src(i + 1, j, k) * Aplusi(i, j, k) +
+                   src(i, j + 1, k) * Aplusj(i, j, k) + src(i, j, k + 1) * Aplusk(i, j, k) +
+                   src(i - 1, j, k) * Aminusi(i, j, k) + src(i, j - 1, k) * Aminusj(i, j, k) +
+                   src(i, j, k - 1) * Aminusk(i, j, k);
+
+    dst(i, j, k) += srcV(i, j + 1, k) * (*matrixA[7])(i, j, k) +
+                    srcV(i - 1, j + 1, k) * (*matrixA[8])(i, j, k) +
+                    srcV(i, j, k) * (*matrixA[9])(i, j, k) +
+                    srcV(i - 1, j, k) * (*matrixA[10])(i, j, k) +
+                    srcW(i, j, k + 1) * (*matrixA[11])(i, j, k) +
+                    srcW(i - 1, j, k + 1) * (*matrixA[12])(i, j, k) +
+                    srcW(i, j, k) * (*matrixA[13])(i, j, k) +
+                    srcW(i - 1, j, k) * (*matrixA[14])(i, j, k);
+  }
+  inline const FlagGrid &getArg0()
+  {
+    return flags;
+  }
+  typedef FlagGrid type0;
+  inline Grid<Real> &getArg1()
+  {
+    return dst;
+  }
+  typedef Grid<Real> type1;
+  inline const Grid<Real> &getArg2()
+  {
+    return src;
+  }
+  typedef Grid<Real> type2;
+  inline const std::vector<Grid<Real> *> &getArg3()
+  {
+    return matrixA;
+  }
+  typedef std::vector<Grid<Real> *> type3;
+  inline const std::vector<Grid<Real> *> &getArg4()
+  {
+    return vecRhs;
+  }
+  typedef std::vector<Grid<Real> *> type4;
   void runMessage()
   {
-    debMsg("Executing kernel ApplyMatrix2D ", 3);
+    debMsg("Executing kernel ApplyMatrixViscosityU ", 3);
     debMsg("Kernel range"
                << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
            4);
   };
   void operator()(const tbb::blocked_range<IndexInt> &__r) const
   {
-    for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
-      op(idx, flags, dst, src, A0, Ai, Aj, Ak);
+    const int _maxX = maxX;
+    const int _maxY = maxY;
+    if (maxZ > 1) {
+      for (int k = __r.begin(); k != (int)__r.end(); k++)
+        for (int j = 1; j < _maxY; j++)
+          for (int i = 1; i < _maxX; i++)
+            op(i, j, k, flags, dst, src, matrixA, vecRhs);
+    }
+    else {
+      const int k = 0;
+      for (int j = __r.begin(); j != (int)__r.end(); j++)
+        for (int i = 1; i < _maxX; i++)
+          op(i, j, k, flags, dst, src, matrixA, vecRhs);
+    }
   }
   void run()
   {
-    tbb::parallel_for(tbb::blocked_range<IndexInt>(0, size), *this);
+    if (maxZ > 1)
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
+    else
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
   }
   const FlagGrid &flags;
   Grid<Real> &dst;
   const Grid<Real> &src;
-  Grid<Real> &A0;
-  Grid<Real> &Ai;
-  Grid<Real> &Aj;
-  Grid<Real> &Ak;
+  const std::vector<Grid<Real> *> matrixA;
+  const std::vector<Grid<Real> *> vecRhs;
 };
 
+struct ApplyMatrixViscosityV : public KernelBase {
+  ApplyMatrixViscosityV(const FlagGrid &flags,
+                        Grid<Real> &dst,
+                        const Grid<Real> &src,
+                        const std::vector<Grid<Real> *> matrixA,
+                        const std::vector<Grid<Real> *> vecRhs)
+      : KernelBase(&flags, 1), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(int i,
+                 int j,
+                 int k,
+                 const FlagGrid &flags,
+                 Grid<Real> &dst,
+                 const Grid<Real> &src,
+                 const std::vector<Grid<Real> *> matrixA,
+                 const std::vector<Grid<Real> *> vecRhs) const
+  {
+    if (matrixA.size() != 15)
+      errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step");
+    Grid<Real> &A0 = *matrixA[0];
+    Grid<Real> &Aplusi = *matrixA[1];
+    Grid<Real> &Aplusj = *matrixA[2];
+    Grid<Real> &Aplusk = *matrixA[3];
+    Grid<Real> &Aminusi = *matrixA[4];
+    Grid<Real> &Aminusj = *matrixA[5];
+    Grid<Real> &Aminusk = *matrixA[6];
+
+    if (vecRhs.size() != 2)
+      errMsg("ConjugatedGrad: Invalid rhs vector in apply matrix step");
+    Grid<Real> &srcU = *vecRhs[0];
+    Grid<Real> &srcW = *vecRhs[1];
+
+    dst(i, j, k) = src(i, j, k) * A0(i, j, k) + src(i + 1, j, k) * Aplusi(i, j, k) +
+                   src(i, j + 1, k) * Aplusj(i, j, k) + src(i, j, k + 1) * Aplusk(i, j, k) +
+                   src(i - 1, j, k) * Aminusi(i, j, k) + src(i, j - 1, k) * Aminusj(i, j, k) +
+                   src(i, j, k - 1) * Aminusk(i, j, k);
+
+    dst(i, j, k) += srcU(i + 1, j, k) * (*matrixA[7])(i, j, k) +
+                    srcU(i + 1, j - 1, k) * (*matrixA[8])(i, j, k) +
+                    srcU(i, j, k) * (*matrixA[9])(i, j, k) +
+                    srcU(i, j - 1, k) * (*matrixA[10])(i, j, k) +
+                    srcW(i, j, k + 1) * (*matrixA[11])(i, j, k) +
+                    srcW(i, j - 1, k + 1) * (*matrixA[12])(i, j, k) +
+                    srcW(i, j, k) * (*matrixA[13])(i, j, k) +
+                    srcW(i, j - 1, k) * (*matrixA[14])(i, j, k);
+  }
+  inline const FlagGrid &getArg0()
+  {
+    return flags;
+  }
+  typedef FlagGrid type0;
+  inline Grid<Real> &getArg1()
+  {
+    return dst;
+  }
+  typedef Grid<Real> type1;
+  inline const Grid<Real> &getArg2()
+  {
+    return src;
+  }
+  typedef Grid<Real> type2;
+  inline const std::vector<Grid<Real> *> &getArg3()
+  {
+    return matrixA;
+  }
+  typedef std::vector<Grid<Real> *> type3;
+  inline const std::vector<Grid<Real> *> &getArg4()
+  {
+    return vecRhs;
+  }
+  typedef std::vector<Grid<Real> *> type4;
+  void runMessage()
+  {
+    debMsg("Executing kernel ApplyMatrixViscosityV ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r) const
+  {
+    const int _maxX = maxX;
+    const int _maxY = maxY;
+    if (maxZ > 1) {
+      for (int k = __r.begin(); k != (int)__r.end(); k++)
+        for (int j = 1; j < _maxY; j++)
+          for (int i = 1; i < _maxX; i++)
+            op(i, j, k, flags, dst, src, matrixA, vecRhs);
+    }
+    else {
+      const int k = 0;
+      for (int j = __r.begin(); j != (int)__r.end(); j++)
+        for (int i = 1; i < _maxX; i++)
+          op(i, j, k, flags, dst, src, matrixA, vecRhs);
+    }
+  }
+  void run()
+  {
+    if (maxZ > 1)
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
+    else
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
+  }
+  const FlagGrid &flags;
+  Grid<Real> &dst;
+  const Grid<Real> &src;
+  const std::vector<Grid<Real> *> matrixA;
+  const std::vector<Grid<Real> *> vecRhs;
+};
+
+struct ApplyMatrixViscosityW : public KernelBase {
+  ApplyMatrixViscosityW(const FlagGrid &flags,
+                        Grid<Real> &dst,
+                        const Grid<Real> &src,
+                        const std::vector<Grid<Real> *> matrixA,
+                        const std::vector<Grid<Real> *> vecRhs)
+      : KernelBase(&flags, 1), flags(flags), dst(dst), src(src), matrixA(matrixA), vecRhs(vecRhs)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(int i,
+                 int j,
+                 int k,
+                 const FlagGrid &flags,
+                 Grid<Real> &dst,
+                 const Grid<Real> &src,
+                 const std::vector<Grid<Real> *> matrixA,
+                 const std::vector<Grid<Real> *> vecRhs) const
+  {
+    if (matrixA.size() != 15)
+      errMsg("ConjugatedGrad: Invalid A matrix in apply matrix step");
+    Grid<Real> &A0 = *matrixA[0];
+    Grid<Real> &Aplusi = *matrixA[1];
+    Grid<Real> &Aplusj = *matrixA[2];
+    Grid<Real> &Aplusk = *matrixA[3];
+    Grid<Real> &Aminusi = *matrixA[4];
+    Grid<Real> &Aminusj = *matrixA[5];
+    Grid<Real> &Aminusk = *matrixA[6];
+
+    if (vecRhs.size() != 2)
+      errMsg("ConjugatedGrad: Invalid rhs vector in apply matrix step");
+    Grid<Real> &srcU = *vecRhs[0];
+    Grid<Real> &srcV = *vecRhs[1];
+
+    dst(i, j, k) = src(i, j, k) * A0(i, j, k) + src(i + 1, j, k) * Aplusi(i, j, k) +
+                   src(i, j + 1, k) * Aplusj(i, j, k) + src(i, j, k + 1) * Aplusk(i, j, k) +
+                   src(i - 1, j, k) * Aminusi(i, j, k) + src(i, j - 1, k) * Aminusj(i, j, k) +
+                   src(i, j, k - 1) * Aminusk(i, j, k);
+
+    dst(i, j, k) += srcU(i + 1, j, k) * (*matrixA[7])(i, j, k) +
+                    srcU(i + 1, j, k - 1) * (*matrixA[8])(i, j, k) +
+                    srcU(i, j, k) * (*matrixA[9])(i, j, k) +
+                    srcU(i, j, k - 1) * (*matrixA[10])(i, j, k) +
+                    srcV(i, j + 1, k) * (*matrixA[11])(i, j, k) +
+                    srcV(i, j + 1, k - 1) * (*matrixA[12])(i, j, k) +
+                    srcV(i, j, k) * (*matrixA[13])(i, j, k) +
+                    srcV(i, j, k - 1) * (*matrixA[14])(i, j, k);
+  }
+  inline const FlagGrid &getArg0()
+  {
+    return flags;
+  }
+  typedef FlagGrid type0;
+  inline Grid<Real> &getArg1()
+  {
+    return dst;
+  }
+  typedef Grid<Real> type1;
+  inline const Grid<Real> &getArg2()
+  {
+    return src;
+  }
+  typedef Grid<Real> type2;
+  inline const std::vector<Grid<Real> *> &getArg3()
+  {
+    return matrixA;
+  }
+  typedef std::vector<Grid<Real> *> type3;
+  inline const std::vector<Grid<Real> *> &getArg4()
+  {
+    return vecRhs;
+  }
+  typedef std::vector<Grid<Real> *> type4;
+  void runMessage()
+  {
+    debMsg("Executing kernel ApplyMatrixViscosityW ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r) const
+  {
+    const int _maxX = maxX;
+    const int _maxY = maxY;
+    if (maxZ > 1) {
+      for (int k = __r.begin(); k != (int)__r.end(); k++)
+        for (int j = 1; j < _maxY; j++)
+          for (int i = 1; i < _maxX; i++)
+            op(i, j, k, flags, dst, src, matrixA, vecRhs);
+    }
+    else {
+      const int k = 0;
+      for (int j = __r.begin(); j != (int)__r.end(); j++)
+        for (int i = 1; i < _maxX; i++)
+          op(i, j, k, flags, dst, src, matrixA, vecRhs);
+    }
+  }
+  void run()
+  {
+    if (maxZ > 1)
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
+    else
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
+  }
+  const FlagGrid &flags;
+  Grid<Real> &dst;
+  const Grid<Real> &src;
+  const std::vector<Grid<Real> *> matrixA;
+  const std::vector<Grid<Real> *> vecRhs;
+};
+
+/* NOTE: Use this template for new matrix application kernels
+
+//! Template for matrix application kernels
+KERNEL()
+void ApplyMatrixTemplate (const FlagGrid& flags, Grid<Real>& dst, const Grid<Real>& src,
+  const std::vector<Grid<Real> *> matrixA, const std::vector<Grid<Real> *> vecRhs)
+{
+  // The kernel must define how to use the grids from the matrixA and vecRhs lists
+}
+
+*/
+
 //! Kernel: Construct the matrix for the poisson equation
 
 struct MakeLaplaceMatrix : public KernelBase {
diff --git a/extern/mantaflow/preprocessed/general.h b/extern/mantaflow/preprocessed/general.h
index 7a840517cef..50eac71e87e 100644
--- a/extern/mantaflow/preprocessed/general.h
+++ b/extern/mantaflow/preprocessed/general.h
@@ -42,7 +42,7 @@ inline void updateQtGui(bool full, int frame, float time, const std::string &cur
 #  ifdef _DEBUG
 #    define DEBUG 1
 #  endif  // _DEBUG
-#endif    // DEBUG
+#endif  // DEBUG
 
 // Standard exception
 class Error : public std::exception {
@@ -242,6 +242,39 @@ inline bool c_isnan(float c)
   return d != d;
 }
 
+//! Swap so that a<b
+template<class T> inline void sort(T &a, T &b)
+{
+  if (a > b)
+    std::swap(a, b);
+}
+
+//! Swap so that a<b<c
+template<class T> inline void sort(T &a, T &b, T &c)
+{
+  if (a > b)
+    std::swap(a, b);
+  if (a > c)
+    std::swap(a, c);
+  if (b > c)
+    std::swap(b, c);
+}
+
+//! Swap so that a<b<c<d
+template<class T> inline void sort(T &a, T &b, T &c, T &d)
+{
+  if (a > b)
+    std::swap(a, b);
+  if (c > d)
+    std::swap(c, d);
+  if (a > c)
+    std::swap(a, c);
+  if (b > d)
+    std::swap(b, d);
+  if (b > c)
+    std::swap(b, c);
+}
+
 }  // namespace Manta
 
 #endif
diff --git a/extern/mantaflow/preprocessed/gitinfo.h b/extern/mantaflow/preprocessed/gitinfo.h
index a0590e6a0b0..afc34a00d3d 100644
--- a/extern/mantaflow/preprocessed/gitinfo.h
+++ b/extern/mantaflow/preprocessed/gitinfo.h
@@ -1,3 +1,3 @@
 
 
-#define MANTA_GIT_VERSION "commit 327917cd59b03bef3a953b5f58fc1637b3a83e01"
+#define MANTA_GIT_VERSION "commit e2285cb9bc492987f728123be6cfc1fe11fe73d6"
diff --git a/extern/mantaflow/preprocessed/plugin/extforces.cpp b/extern/mantaflow/preprocessed/plugin/extforces.cpp
index 558008afe3b..88935fa7ae9 100644
--- a/extern/mantaflow/preprocessed/plugin/extforces.cpp
+++ b/extern/mantaflow/preprocessed/plugin/extforces.cpp
@@ -1135,26 +1135,27 @@ struct KnAddForceIfLower : public KernelBase {
     if (!curFluid && !curEmpty)
       return;
 
+    Real minVal, maxVal, sum;
     if (flags.isFluid(i - 1, j, k) || (curFluid && flags.isEmpty(i - 1, j, k))) {
       Real forceMACX = 0.5 * (force(i - 1, j, k).x + force(i, j, k).x);
-      Real min = std::min(vel(i, j, k).x, forceMACX);
-      Real max = std::max(vel(i, j, k).x, forceMACX);
-      Real sum = vel(i, j, k).x + forceMACX;
-      vel(i, j, k).x = (forceMACX > 0) ? std::min(sum, max) : std::max(sum, min);
+      minVal = min(vel(i, j, k).x, forceMACX);
+      maxVal = max(vel(i, j, k).x, forceMACX);
+      sum = vel(i, j, k).x + forceMACX;
+      vel(i, j, k).x = (forceMACX > 0) ? min(sum, maxVal) : max(sum, minVal);
     }
     if (flags.isFluid(i, j - 1, k) || (curFluid && flags.isEmpty(i, j - 1, k))) {
       Real forceMACY = 0.5 * (force(i, j - 1, k).y + force(i, j, k).y);
-      Real min = std::min(vel(i, j, k).y, forceMACY);
-      Real max = std::max(vel(i, j, k).y, forceMACY);
-      Real sum = vel(i, j, k).y + forceMACY;
-      vel(i, j, k).y = (forceMACY > 0) ? std::min(sum, max) : std::max(sum, min);
+      minVal = min(vel(i, j, k).y, forceMACY);
+      maxVal = max(vel(i, j, k).y, forceMACY);
+      sum = vel(i, j, k).y + forceMACY;
+      vel(i, j, k).y = (forceMACY > 0) ? min(sum, maxVal) : max(sum, minVal);
     }
     if (vel.is3D() && (flags.isFluid(i, j, k - 1) || (curFluid && flags.isEmpty(i, j, k - 1)))) {
       Real forceMACZ = 0.5 * (force(i, j, k - 1).z + force(i, j, k).z);
-      Real min = std::min(vel(i, j, k).z, forceMACZ);
-      Real max = std::max(vel(i, j, k).z, forceMACZ);
-      Real sum = vel(i, j, k).z + forceMACZ;
-      vel(i, j, k).z = (forceMACZ > 0) ? std::min(sum, max) : std::max(sum, min);
+      minVal = min(vel(i, j, k).z, forceMACZ);
+      maxVal = max(vel(i, j, k).z, forceMACZ);
+      sum = vel(i, j, k).z + forceMACZ;
+      vel(i, j, k).z = (forceMACZ > 0) ? min(sum, maxVal) : max(sum, minVal);
     }
   }
   inline const FlagGrid &getArg0()
diff --git a/extern/mantaflow/preprocessed/plugin/pressure.cpp b/extern/mantaflow/preprocessed/plugin/pressure.cpp
index 1100a58db47..593aeb16859 100644
--- a/extern/mantaflow/preprocessed/plugin/pressure.cpp
+++ b/extern/mantaflow/preprocessed/plugin/pressure.cpp
@@ -1138,11 +1138,15 @@ void solvePressureSystem(Grid<Real> &rhs,
   // note: the last factor increases the max iterations for 2d, which right now can't use a
   // preconditioner
   GridCgInterface *gcg;
-  if (vel.is3D())
-    gcg = new GridCg<ApplyMatrix>(pressure, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
-  else
-    gcg = new GridCg<ApplyMatrix2D>(
-        pressure, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+  vector<Grid<Real> *> matA{&A0, &Ai, &Aj};
+
+  if (vel.is3D()) {
+    matA.push_back(&Ak);
+    gcg = new GridCg<ApplyMatrix>(pressure, rhs, residual, search, flags, tmp, matA);
+  }
+  else {
+    gcg = new GridCg<ApplyMatrix2D>(pressure, rhs, residual, search, flags, tmp, matA);
+  }
 
   gcg->setAccuracy(cgAccuracy);
   gcg->setUseL2Norm(useL2Norm);
diff --git a/extern/mantaflow/preprocessed/plugin/viscosity.cpp b/extern/mantaflow/preprocessed/plugin/viscosity.cpp
new file mode 100644
index 00000000000..5dc478b13dd
--- /dev/null
+++ b/extern/mantaflow/preprocessed/plugin/viscosity.cpp
@@ -0,0 +1,1430 @@
+
+
+// DO NOT EDIT !
+// This file is generated using the MantaFlow preprocessor (prep generate).
+
+/******************************************************************************
+ *
+ * MantaFlow fluid solver framework
+ * Copyright 2020 Sebastian Barschkis, Nils Thuerey
+ *
+ * This program is free software, distributed under the terms of the
+ * Apache License, Version 2.0
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Accurate Viscous Free Surfaces for Buckling, Coiling, and Rotating Liquids
+ * Batty et al., SCA 2008
+ *
+ ******************************************************************************/
+
+#include "conjugategrad.h"
+#include "general.h"
+#include "grid.h"
+#include "vectorbase.h"
+
+#include <chrono>
+
+#if OPENMP == 1 || TBB == 1
+#  define ENABLE_PARALLEL 0
+#endif
+
+#if ENABLE_PARALLEL == 1
+#  include <thread>
+#  include <algorithm>
+
+static const int manta_num_threads = std::thread::hardware_concurrency();
+
+#  define parallel_block \
+    { \
+      std::vector<std::thread> threads; \
+      {
+
+#  define do_parallel threads.push_back( std::thread([&]() {
+#  define do_end \
+    } ) );
+
+#  define block_end \
+    } \
+    for (auto &thread : threads) { \
+      thread.join(); \
+    } \
+    }
+
+#endif
+
+#define FOR_INT_IJK(num) \
+  for (int k_off = 0; k_off < num; ++k_off) \
+    for (int j_off = 0; j_off < num; ++j_off) \
+      for (int i_off = 0; i_off < num; ++i_off)
+
+using namespace std;
+
+namespace Manta {
+
+//! Assumes phi0<0 and phi1>=0, phi2>=0, and phi3>=0 or vice versa.
+//! In particular, phi0 must not equal any of phi1, phi2 or phi3.
+static Real sortedTetFraction(Real phi0, Real phi1, Real phi2, Real phi3)
+{
+  return phi0 * phi0 * phi0 / ((phi0 - phi1) * (phi0 - phi2) * (phi0 - phi3));
+}
+
+//! Assumes phi0<0, phi1<0, and phi2>=0, and phi3>=0 or vice versa.
+//! In particular, phi0 and phi1 must not equal any of phi2 and phi3.
+static Real sortedPrismFraction(Real phi0, Real phi1, Real phi2, Real phi3)
+{
+  Real a = phi0 / (phi0 - phi2);
+  Real b = phi0 / (phi0 - phi3);
+  Real c = phi1 / (phi1 - phi3);
+  Real d = phi1 / (phi1 - phi2);
+  return a * b * (1 - d) + b * (1 - c) * d + c * d;
+}
+
+Real volumeFraction(Real phi0, Real phi1, Real phi2, Real phi3)
+{
+  sort(phi0, phi1, phi2, phi3);
+  if (phi3 <= 0)
+    return 1;
+  else if (phi2 <= 0)
+    return 1 - sortedTetFraction(phi3, phi2, phi1, phi0);
+  else if (phi1 <= 0)
+    return sortedPrismFraction(phi0, phi1, phi2, phi3);
+  else if (phi0 <= 0)
+    return sortedTetFraction(phi0, phi1, phi2, phi3);
+  else
+    return 0;
+}
+
+//! The average of the two possible decompositions of the cube into five tetrahedra.
+Real volumeFraction(Real phi000,
+                    Real phi100,
+                    Real phi010,
+                    Real phi110,
+                    Real phi001,
+                    Real phi101,
+                    Real phi011,
+                    Real phi111)
+{
+  return (volumeFraction(phi000, phi001, phi101, phi011) +
+          volumeFraction(phi000, phi101, phi100, phi110) +
+          volumeFraction(phi000, phi010, phi011, phi110) +
+          volumeFraction(phi101, phi011, phi111, phi110) +
+          2 * volumeFraction(phi000, phi011, phi101, phi110) +
+          volumeFraction(phi100, phi101, phi001, phi111) +
+          volumeFraction(phi100, phi001, phi000, phi010) +
+          volumeFraction(phi100, phi110, phi111, phi010) +
+          volumeFraction(phi001, phi111, phi011, phi010) +
+          2 * volumeFraction(phi100, phi111, phi001, phi010)) /
+         12;
+}
+
+//! Kernel loop over grid with 2x base resolution!
+
+struct KnEstimateVolumeFraction : public KernelBase {
+  KnEstimateVolumeFraction(Grid<Real> &volumes,
+                           const Grid<Real> &phi,
+                           const Vec3 &startCentre,
+                           const Real dx)
+      : KernelBase(&volumes, 0), volumes(volumes), phi(phi), startCentre(startCentre), dx(dx)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(int i,
+                 int j,
+                 int k,
+                 Grid<Real> &volumes,
+                 const Grid<Real> &phi,
+                 const Vec3 &startCentre,
+                 const Real dx) const
+  {
+    const Vec3 centre = startCentre + Vec3(i, j, k) * 0.5;
+    const Real offset = 0.5 * dx;
+    const int order = 2;
+
+    Real phi000 = phi.getInterpolatedHi(centre + Vec3(-offset, -offset, -offset), order);
+    Real phi001 = phi.getInterpolatedHi(centre + Vec3(-offset, -offset, +offset), order);
+    Real phi010 = phi.getInterpolatedHi(centre + Vec3(-offset, +offset, -offset), order);
+    Real phi011 = phi.getInterpolatedHi(centre + Vec3(-offset, +offset, +offset), order);
+    Real phi100 = phi.getInterpolatedHi(centre + Vec3(+offset, -offset, -offset), order);
+    Real phi101 = phi.getInterpolatedHi(centre + Vec3(+offset, -offset, +offset), order);
+    Real phi110 = phi.getInterpolatedHi(centre + Vec3(+offset, +offset, -offset), order);
+    Real phi111 = phi.getInterpolatedHi(centre + Vec3(+offset, +offset, +offset), order);
+
+    volumes(i, j, k) = volumeFraction(
+        phi000, phi100, phi010, phi110, phi001, phi101, phi011, phi111);
+  }
+  inline Grid<Real> &getArg0()
+  {
+    return volumes;
+  }
+  typedef Grid<Real> type0;
+  inline const Grid<Real> &getArg1()
+  {
+    return phi;
+  }
+  typedef Grid<Real> type1;
+  inline const Vec3 &getArg2()
+  {
+    return startCentre;
+  }
+  typedef Vec3 type2;
+  inline const Real &getArg3()
+  {
+    return dx;
+  }
+  typedef Real type3;
+  void runMessage()
+  {
+    debMsg("Executing kernel KnEstimateVolumeFraction ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r) const
+  {
+    const int _maxX = maxX;
+    const int _maxY = maxY;
+    if (maxZ > 1) {
+      for (int k = __r.begin(); k != (int)__r.end(); k++)
+        for (int j = 0; j < _maxY; j++)
+          for (int i = 0; i < _maxX; i++)
+            op(i, j, k, volumes, phi, startCentre, dx);
+    }
+    else {
+      const int k = 0;
+      for (int j = __r.begin(); j != (int)__r.end(); j++)
+        for (int i = 0; i < _maxX; i++)
+          op(i, j, k, volumes, phi, startCentre, dx);
+    }
+  }
+  void run()
+  {
+    if (maxZ > 1)
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
+    else
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(0, maxY), *this);
+  }
+  Grid<Real> &volumes;
+  const Grid<Real> &phi;
+  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();