Physics: Improve liquid simulation performance

Move various computations out of a loop.

This can improve performance up to 1.25x - 1.5x depending on the scene.

Pull Request: https://projects.blender.org/blender/blender/pulls/137733

extern/mantaflow/README.blender

@@ -6,3 +6,5 @@ Upstream version: 0.13
 Local modifications:
 * ./patches/local_namespace.diff to support loading MANTA variables into an isolated __main__ name-space.
 * ./patches/fix-computation-errors.patch fix computation errors for normalization functions for 3d vectors by using std::hypot.
+* ./patches/liquid-mesh-performance.patch improve liquid mesh generation by puting calculation of inverse radius outside for loops.
+* ./patches/liquid-performance.patch improve liquid generation (without mesh) by precalculate sum of vectors and put it outside for loop.

extern/mantaflow/patches/liquid-mesh-performance.patch

@@ -0,0 +1,43 @@
+diff --git a/extern/mantaflow/preprocessed/plugin/flip.cpp b/extern/mantaflow/preprocessed/plugin/flip.cpp
+index 8757958d4b0..f5d7147c34f 100644
+--- a/extern/mantaflow/preprocessed/plugin/flip.cpp
++++ b/extern/mantaflow/preprocessed/plugin/flip.cpp
+@@ -1012,6 +1012,7 @@ struct ComputeAveragedLevelsetWeight : public KernelBase {
+                  const ParticleIndexSystem &indexSys,
+                  LevelsetGrid &phi,
+                  const Real radius,
++                 const Real sradiusInv,
+                  const ParticleDataImpl<int> *ptype,
+                  const int exclude,
+                  Grid<Vec3> *save_pAcc = nullptr,
+@@ -1020,8 +1021,6 @@ struct ComputeAveragedLevelsetWeight : public KernelBase {
+     const Vec3 gridPos = Vec3(i, j, k) + Vec3(0.5);  // shifted by half cell
+     Real phiv = radius * 1.0;                        // outside
+ 
+-    // loop over neighborhood, similar to ComputeUnionLevelsetPindex
+-    const Real sradiusInv = 1. / (4. * radius * radius);
+     const int r = int(radius) + 1;
+     // accumulators
+     Real wacc = 0.;
+@@ -1120,17 +1119,19 @@ struct ComputeAveragedLevelsetWeight : public KernelBase {
+   {
+     const int _maxX = maxX;
+     const int _maxY = maxY;
++    // loop over neighborhood, similar to ComputeUnionLevelsetPindex
++    const Real sradiusInv = 1. / (4. * radius * radius);
+     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, parts, index, indexSys, phi, radius, ptype, exclude, save_pAcc, save_rAcc);
++            op(i, j, k, parts, index, indexSys, phi, radius, sradiusInv, ptype, exclude, save_pAcc, save_rAcc);
+     }
+     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, parts, index, indexSys, phi, radius, ptype, exclude, save_pAcc, save_rAcc);
++          op(i, j, k, parts, index, indexSys, phi, radius, sradiusInv, ptype, exclude, save_pAcc, save_rAcc);
+     }
+   }
+   void run()

extern/mantaflow/patches/liquid-performance.patch

@@ -0,0 +1,73 @@
+diff --git a/extern/mantaflow/preprocessed/fastmarch.cpp b/extern/mantaflow/preprocessed/fastmarch.cpp
+index 31e43483b49..9856d84a8b1 100644
+--- a/extern/mantaflow/preprocessed/fastmarch.cpp
++++ b/extern/mantaflow/preprocessed/fastmarch.cpp
+@@ -306,25 +306,24 @@ struct knExtrapolateMACSimple : public KernelBase {
+                  const int d,
+                  const int c) const
+   {
+-    static const Vec3i nb[6] = {Vec3i(1, 0, 0),
+-                                Vec3i(-1, 0, 0),
+-                                Vec3i(0, 1, 0),
+-                                Vec3i(0, -1, 0),
+-                                Vec3i(0, 0, 1),
+-                                Vec3i(0, 0, -1)};
+-    const int dim = (vel.is3D() ? 3 : 2);
+-
+     if (tmp(i, j, k) != 0)
+       return;
++    static const Vec3i nb[6] = {Vec3i(i+1, j, k),
++                                Vec3i(i-1, j, k),
++                                Vec3i(i, j+1, k),
++                                Vec3i(i, j-1, k),
++                                Vec3i(i, j, k+1),
++                                Vec3i(i, j, k-1)};
++    const int dim = (vel.is3D() ? 3 : 2);
+ 
+     // copy from initialized neighbors
+     Vec3i p(i, j, k);
+     int nbs = 0;
+     Real avgVel = 0.;
+     for (int n = 0; n < 2 * dim; ++n) {
+-      if (tmp(p + nb[n]) == d) {
++      if (tmp(nb[n]) == d) {
+         // vel(p)[c] = (c+1.)*0.1;
+-        avgVel += vel(p + nb[n])[c];
++        avgVel += vel(nb[n])[c];
+         nbs++;
+       }
+     }
+@@ -714,24 +713,23 @@ struct knExtrapolateMACFromWeight : public KernelBase {
+                  const int d,
+                  const int c) const
+   {
+-    static const Vec3i nb[6] = {Vec3i(1, 0, 0),
+-                                Vec3i(-1, 0, 0),
+-                                Vec3i(0, 1, 0),
+-                                Vec3i(0, -1, 0),
+-                                Vec3i(0, 0, 1),
+-                                Vec3i(0, 0, -1)};
+-    const int dim = (vel.is3D() ? 3 : 2);
+-
+     if (weight(i, j, k)[c] != 0)
+       return;
++    static const Vec3i nb[6] = {Vec3i(i+1, j, k),
++                                Vec3i(i-1, j, k),
++                                Vec3i(i, j+1, k),
++                                Vec3i(i, j-1, k),
++                                Vec3i(i, j, k+1),
++                                Vec3i(i, j, k-1)};
++    const int dim = (vel.is3D() ? 3 : 2);
+ 
+     // copy from initialized neighbors
+     Vec3i p(i, j, k);
+     int nbs = 0;
+     Real avgVel = 0.;
+     for (int n = 0; n < 2 * dim; ++n) {
+-      if (weight(p + nb[n])[c] == d) {
+-        avgVel += vel(p + nb[n])[c];
++      if (weight(nb[n])[c] == d) {
++        avgVel += vel(nb[n])[c];
+         nbs++;
+       }
+     }

extern/mantaflow/preprocessed/fastmarch.cpp

@@ -306,25 +306,24 @@ struct knExtrapolateMACSimple : public KernelBase {
                  const int d,
                  const int c) const
   {
-    static const Vec3i nb[6] = {Vec3i(1, 0, 0),
-                                Vec3i(-1, 0, 0),
-                                Vec3i(0, 1, 0),
-                                Vec3i(0, -1, 0),
-                                Vec3i(0, 0, 1),
-                                Vec3i(0, 0, -1)};
-    const int dim = (vel.is3D() ? 3 : 2);
-
     if (tmp(i, j, k) != 0)
       return;
+    static const Vec3i nb[6] = {Vec3i(i+1, j, k),
+                                Vec3i(i-1, j, k),
+                                Vec3i(i, j+1, k),
+                                Vec3i(i, j-1, k),
+                                Vec3i(i, j, k+1),
+                                Vec3i(i, j, k-1)};
+    const int dim = (vel.is3D() ? 3 : 2);
 
     // copy from initialized neighbors
     Vec3i p(i, j, k);
     int nbs = 0;
     Real avgVel = 0.;
     for (int n = 0; n < 2 * dim; ++n) {
-      if (tmp(p + nb[n]) == d) {
+      if (tmp(nb[n]) == d) {
         // vel(p)[c] = (c+1.)*0.1;
-        avgVel += vel(p + nb[n])[c];
+        avgVel += vel(nb[n])[c];
         nbs++;
       }
     }
@@ -714,24 +713,23 @@ struct knExtrapolateMACFromWeight : public KernelBase {
                  const int d,
                  const int c) const
   {
-    static const Vec3i nb[6] = {Vec3i(1, 0, 0),
-                                Vec3i(-1, 0, 0),
-                                Vec3i(0, 1, 0),
-                                Vec3i(0, -1, 0),
-                                Vec3i(0, 0, 1),
-                                Vec3i(0, 0, -1)};
-    const int dim = (vel.is3D() ? 3 : 2);
-
     if (weight(i, j, k)[c] != 0)
       return;
+    static const Vec3i nb[6] = {Vec3i(i+1, j, k),
+                                Vec3i(i-1, j, k),
+                                Vec3i(i, j+1, k),
+                                Vec3i(i, j-1, k),
+                                Vec3i(i, j, k+1),
+                                Vec3i(i, j, k-1)};
+    const int dim = (vel.is3D() ? 3 : 2);
 
     // copy from initialized neighbors
     Vec3i p(i, j, k);
     int nbs = 0;
     Real avgVel = 0.;
     for (int n = 0; n < 2 * dim; ++n) {
-      if (weight(p + nb[n])[c] == d) {
+      if (weight(nb[n])[c] == d) {
-        avgVel += vel(p + nb[n])[c];
+        avgVel += vel(nb[n])[c];
         nbs++;
       }
     }

extern/mantaflow/preprocessed/plugin/flip.cpp

@@ -1012,6 +1012,7 @@ struct ComputeAveragedLevelsetWeight : public KernelBase {
                  const ParticleIndexSystem &indexSys,
                  LevelsetGrid &phi,
                  const Real radius,
+                 const Real sradiusInv,
                  const ParticleDataImpl<int> *ptype,
                  const int exclude,
                  Grid<Vec3> *save_pAcc = nullptr,
@@ -1021,7 +1022,6 @@ struct ComputeAveragedLevelsetWeight : public KernelBase {
     Real phiv = radius * 1.0;                        // outside
 
     // loop over neighborhood, similar to ComputeUnionLevelsetPindex
-    const Real sradiusInv = 1. / (4. * radius * radius);
     const int r = int(radius) + 1;
     // accumulators
     Real wacc = 0.;
@@ -1120,17 +1120,18 @@ struct ComputeAveragedLevelsetWeight : public KernelBase {
   {
     const int _maxX = maxX;
     const int _maxY = maxY;
+    const Real sradiusInv = 1. / (4. * radius * radius);
     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, parts, index, indexSys, phi, radius, ptype, exclude, save_pAcc, save_rAcc);
+            op(i, j, k, parts, index, indexSys, phi, radius, sradiusInv, ptype, exclude, save_pAcc, save_rAcc);
     }
     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, parts, index, indexSys, phi, radius, ptype, exclude, save_pAcc, save_rAcc);
+          op(i, j, k, parts, index, indexSys, phi, radius, sradiusInv, ptype, exclude, save_pAcc, save_rAcc);
     }
   }
   void run()