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.

4 files changed, 53 insertions, 4 deletions

diff --git a/intern/mantaflow/intern/MANTA_main.cpp b/intern/mantaflow/intern/MANTA_main.cpp
index d49915e2452..d5d41e71f22 100644
--- a/intern/mantaflow/intern/MANTA_main.cpp
+++ b/intern/mantaflow/intern/MANTA_main.cpp
@@ -72,6 +72,7 @@ MANTA::MANTA(int *res, FluidModifierData *fmd)
   mUsingFractions = (fds->flags & FLUID_DOMAIN_USE_FRACTIONS) && mUsingLiquid;
   mUsingMesh = (fds->flags & FLUID_DOMAIN_USE_MESH) && mUsingLiquid;
   mUsingDiffusion = (fds->flags & FLUID_DOMAIN_USE_DIFFUSION) && mUsingLiquid;
+  mUsingViscosity = (fds->flags & FLUID_DOMAIN_USE_VISCOSITY) && mUsingLiquid;
   mUsingMVel = (fds->flags & FLUID_DOMAIN_USE_SPEED_VECTORS) && mUsingLiquid;
   mUsingGuiding = (fds->flags & FLUID_DOMAIN_USE_GUIDE);
   mUsingDrops = (fds->particle_type & FLUID_DOMAIN_PARTICLE_SPRAY) && mUsingLiquid;
@@ -221,6 +222,10 @@ MANTA::MANTA(int *res, FluidModifierData *fmd)
       initSuccess &= initLiquidMesh();
     }
 
+    if (mUsingViscosity) {
+      initSuccess &= initLiquidViscosity();
+    }
+
     if (mUsingDiffusion) {
       initSuccess &= initCurvature();
     }
@@ -440,6 +445,17 @@ bool MANTA::initLiquidMesh(FluidModifierData *fmd)
   return runPythonString(pythonCommands);
 }
 
+bool MANTA::initLiquidViscosity(FluidModifierData *fmd)
+{
+  vector<string> pythonCommands;
+  string tmpString = fluid_variables_viscosity + fluid_solver_viscosity + liquid_alloc_viscosity;
+  string finalString = parseScript(tmpString, fmd);
+  pythonCommands.push_back(finalString);
+
+  mUsingViscosity = true;
+  return runPythonString(pythonCommands);
+}
+
 bool MANTA::initCurvature(FluidModifierData *fmd)
 {
   std::vector<std::string> pythonCommands;
@@ -871,8 +887,10 @@ void MANTA::initializeRNAMap(FluidModifierData *fmd)
   mRNAMap["CACHE_DIR"] = cacheDirectory;
   mRNAMap["COMPRESSION_OPENVDB"] = vdbCompressionMethod;
   mRNAMap["PRECISION_OPENVDB"] = vdbPrecisionHalf;
-	mRNAMap["CLIP_OPENVDB"] = to_string(fds->clipping);
+  mRNAMap["CLIP_OPENVDB"] = to_string(fds->clipping);
   mRNAMap["PP_PARTICLE_MAXIMUM"] = to_string(fds->sys_particle_maximum);
+  mRNAMap["USING_VISCOSITY"] = getBooleanString(fds->flags & FLUID_DOMAIN_USE_VISCOSITY);
+  mRNAMap["VISCOSITY_VALUE"] = to_string(fds->viscosity_value);
 
   /* Fluid object names. */
   mRNAMap["NAME_FLAGS"] = FLUID_NAME_FLAGS;
@@ -1728,6 +1746,7 @@ bool MANTA::exportLiquidScript(FluidModifierData *fmd)
   bool guiding = fds->active_fields & FLUID_DOMAIN_ACTIVE_GUIDE;
   bool invel = fds->active_fields & FLUID_DOMAIN_ACTIVE_INVEL;
   bool outflow = fds->active_fields & FLUID_DOMAIN_ACTIVE_OUTFLOW;
+  bool viscosity = fds->flags & FLUID_DOMAIN_USE_VISCOSITY;
 
   string manta_script;
 
@@ -1742,6 +1761,8 @@ bool MANTA::exportLiquidScript(FluidModifierData *fmd)
     manta_script += fluid_variables_particles + liquid_variables_particles;
   if (guiding)
     manta_script += fluid_variables_guiding;
+  if (viscosity)
+    manta_script += fluid_variables_viscosity;
 
   /* Solvers. */
   manta_script += header_solvers + fluid_solver;
@@ -1751,6 +1772,8 @@ bool MANTA::exportLiquidScript(FluidModifierData *fmd)
     manta_script += fluid_solver_particles;
   if (guiding)
     manta_script += fluid_solver_guiding;
+  if (viscosity)
+    manta_script += fluid_solver_viscosity;
 
   /* Grids. */
   manta_script += header_grids + fluid_alloc + liquid_alloc;
@@ -1768,6 +1791,8 @@ bool MANTA::exportLiquidScript(FluidModifierData *fmd)
     manta_script += fluid_alloc_invel;
   if (outflow)
     manta_script += fluid_alloc_outflow;
+  if (viscosity)
+    manta_script += liquid_alloc_viscosity;
 
   /* Domain init. */
   manta_script += header_gridinit + liquid_init_phi;
diff --git a/intern/mantaflow/intern/MANTA_main.h b/intern/mantaflow/intern/MANTA_main.h
index 163b168e43d..7129a545243 100644
--- a/intern/mantaflow/intern/MANTA_main.h
+++ b/intern/mantaflow/intern/MANTA_main.h
@@ -67,6 +67,7 @@ struct MANTA {
   bool initColorsHigh(struct FluidModifierData *fmd = nullptr);
   bool initLiquid(FluidModifierData *fmd = nullptr);
   bool initLiquidMesh(FluidModifierData *fmd = nullptr);
+  bool initLiquidViscosity(FluidModifierData *fmd = nullptr);
   bool initObstacle(FluidModifierData *fmd = nullptr);
   bool initCurvature(FluidModifierData *fmd = nullptr);
   bool initGuiding(FluidModifierData *fmd = nullptr);
@@ -757,6 +758,7 @@ struct MANTA {
   bool mUsingNoise;
   bool mUsingMesh;
   bool mUsingDiffusion;
+  bool mUsingViscosity;
   bool mUsingMVel;
   bool mUsingLiquid;
   bool mUsingSmoke;
diff --git a/intern/mantaflow/intern/strings/fluid_script.h b/intern/mantaflow/intern/strings/fluid_script.h
index 8c9025dd435..73e8552d260 100644
--- a/intern/mantaflow/intern/strings/fluid_script.h
+++ b/intern/mantaflow/intern/strings/fluid_script.h
@@ -79,6 +79,11 @@ const std::string fluid_solver_guiding =
 mantaMsg('Solver guiding')\n\
 sg$ID$ = Solver(name='solver_guiding$ID$', gridSize=gs_sg$ID$)\n";
 
+const std::string fluid_solver_viscosity =
+    "\n\
+mantaMsg('Solver viscosity')\n\
+sv$ID$ = Solver(name='solver_viscosity$ID$', gridSize=gs_sv$ID$, dim=dim_s$ID$)\n";
+
 //////////////////////////////////////////////////////////////////////
 // VARIABLES
 //////////////////////////////////////////////////////////////////////
@@ -133,7 +138,7 @@ end_frame_s$ID$     = $END_FRAME$\n\
 \n\
 # Fluid diffusion / viscosity\n\
 domainSize_s$ID$ = $FLUID_DOMAIN_SIZE$ # longest domain side in meters\n\
-viscosity_s$ID$ = $FLUID_VISCOSITY$ / (domainSize_s$ID$*domainSize_s$ID$) # kinematic viscosity in m^2/s\n\
+kinViscosity_s$ID$ = $FLUID_VISCOSITY$ / (domainSize_s$ID$*domainSize_s$ID$) # kinematic viscosity in m^2/s\n\
 \n\
 # Factors to convert Blender units to Manta units\n\
 ratioMetersToRes_s$ID$ = float(domainSize_s$ID$) / float(res_s$ID$) # [meters / cells]\n\
@@ -199,6 +204,10 @@ tau_sg$ID$   = 1.0\n\
 sigma_sg$ID$ = 0.99/tau_sg$ID$\n\
 theta_sg$ID$ = 1.0\n";
 
+const std::string fluid_variables_viscosity =
+    "\n\
+gs_sv$ID$    = vec3($RESX$*2, $RESY$*2, $RESZ$*2)\n";
+
 const std::string fluid_with_obstacle =
     "\n\
 using_obstacle_s$ID$ = True\n";
diff --git a/intern/mantaflow/intern/strings/liquid_script.h b/intern/mantaflow/intern/strings/liquid_script.h
index 8e0d113d680..c44727bd47e 100644
--- a/intern/mantaflow/intern/strings/liquid_script.h
+++ b/intern/mantaflow/intern/strings/liquid_script.h
@@ -40,6 +40,8 @@ radiusFactor_s$ID$     = $PARTICLE_RADIUS$\n\
 using_mesh_s$ID$       = $USING_MESH$\n\
 using_final_mesh_s$ID$ = $USING_IMPROVED_MESH$\n\
 using_fractions_s$ID$  = $USING_FRACTIONS$\n\
+using_apic_s$ID$       = $USING_APIC$\n\
+using_viscosity_s$ID$  = $USING_VISCOSITY$\n\
 fracThreshold_s$ID$    = $FRACTIONS_THRESHOLD$\n\
 fracDistance_s$ID$     = $FRACTIONS_DISTANCE$\n\
 flipRatio_s$ID$        = $FLIP_RATIO$\n\
@@ -51,7 +53,7 @@ smoothenNeg_s$ID$      = $MESH_SMOOTHEN_NEG$\n\
 randomness_s$ID$       = $PARTICLE_RANDOMNESS$\n\
 surfaceTension_s$ID$   = $LIQUID_SURFACE_TENSION$\n\
 maxSysParticles_s$ID$  = $PP_PARTICLE_MAXIMUM$\n\
-using_apic_s$ID$       = $USING_APIC$\n";
+viscosityValue_s$ID$   = $VISCOSITY_VALUE$\n";
 
 const std::string liquid_variables_particles =
     "\n\
@@ -135,6 +137,13 @@ liquid_mesh_dict_s$ID$ = { 'lMesh' : mesh_sm$ID$ }\n\
 if using_speedvectors_s$ID$:\n\
     liquid_meshvel_dict_s$ID$ = { 'lVelMesh' : mVel_mesh$ID$ }\n";
 
+const std::string liquid_alloc_viscosity =
+    "\n\
+# Viscosity grids\n\
+volumes_s$ID$   = sv$ID$.create(RealGrid)\n\
+viscosity_s$ID$ = s$ID$.create(RealGrid)\n\
+viscosity_s$ID$.setConst(viscosityValue_s$ID$)\n";
+
 const std::string liquid_alloc_curvature =
     "\n\
 mantaMsg('Liquid alloc curvature')\n\
@@ -306,7 +315,7 @@ def liquid_step_$ID$():\n\
     if using_diffusion_s$ID$:\n\
         mantaMsg('Viscosity')\n\
         # diffusion param for solve = const * dt / dx^2\n\
-        alphaV = viscosity_s$ID$ * s$ID$.timestep * float(res_s$ID$*res_s$ID$)\n\
+        alphaV = kinViscosity_s$ID$ * s$ID$.timestep * float(res_s$ID$*res_s$ID$)\n\
         setWallBcs(flags=flags_s$ID$, vel=vel_s$ID$, obvel=None if using_fractions_s$ID$ else obvel_s$ID$, phiObs=phiObs_s$ID$, fractions=fractions_s$ID$)\n\
         cgSolveDiffusion(flags_s$ID$, vel_s$ID$, alphaV)\n\
         \n\
@@ -315,7 +324,11 @@ def liquid_step_$ID$():\n\
         curvature_s$ID$.clamp(-1.0, 1.0)\n\
     \n\
     setWallBcs(flags=flags_s$ID$, vel=vel_s$ID$, obvel=None if using_fractions_s$ID$ else obvel_s$ID$, phiObs=phiObs_s$ID$, fractions=fractions_s$ID$)\n\
+    if using_viscosity_s$ID$:\n\
+        viscosity_s$ID$.setConst(viscosityValue_s$ID$)\n\
+        applyViscosity(flags=flags_s$ID$, phi=phi_s$ID$, vel=vel_s$ID$, volumes=volumes_s$ID$, viscosity=viscosity_s$ID$)\n\
     \n\
+    setWallBcs(flags=flags_s$ID$, vel=vel_s$ID$, obvel=None if using_fractions_s$ID$ else obvel_s$ID$, phiObs=phiObs_s$ID$, fractions=fractions_s$ID$)\n\
     if using_guiding_s$ID$:\n\
         mantaMsg('Guiding and pressure')\n\
         PD_fluid_guiding(vel=vel_s$ID$, velT=velT_s$ID$, flags=flags_s$ID$, phi=phi_s$ID$, curv=curvature_s$ID$, surfTens=surfaceTension_s$ID$, fractions=fractions_s$ID$, weight=weightGuide_s$ID$, blurRadius=beta_sg$ID$, pressure=pressure_s$ID$, tau=tau_sg$ID$, sigma=sigma_sg$ID$, theta=theta_sg$ID$, zeroPressureFixing=domainClosed_s$ID$)\n\