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.

1 files changed, 1430 insertions, 0 deletions

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