Fluid: Added APIC simulation method

Basic support for velocity updates with the APIC method.

This commit adds APIC to the already existing dropdown menu for the simulation method. The APIC plugin within Mantaflow has been updated to the latest version.

6 files changed, 280 insertions, 155 deletions

diff --git a/extern/mantaflow/preprocessed/plugin/apic.cpp b/extern/mantaflow/preprocessed/plugin/apic.cpp
index 6ff893014c9..100b26a4091 100644
--- a/extern/mantaflow/preprocessed/plugin/apic.cpp
+++ b/extern/mantaflow/preprocessed/plugin/apic.cpp
@@ -6,7 +6,7 @@
 // ----------------------------------------------------------------------------
 //
 // MantaFlow fluid solver framework
-// Copyright 2016-2017 Kiwon Um, Nils Thuerey
+// Copyright 2016-2020 Kiwon Um, Nils Thuerey
 //
 // This program is free software, distributed under the terms of the
 // Apache License, Version 2.0
@@ -21,6 +21,54 @@
 
 namespace Manta {
 
+#define FOR_INT_IJK(num) \
+  for (int i = 0; i < num; ++i) \
+    for (int j = 0; j < num; ++j) \
+      for (int k = 0; k < num; ++k)
+
+static inline IndexInt indexUFace(const Vec3 &pos, const MACGrid &ref)
+{
+  const Vec3i f = toVec3i(pos), c = toVec3i(pos - 0.5);
+  const IndexInt index = f.x * ref.getStrideX() + c.y * ref.getStrideY() + c.z * ref.getStrideZ();
+  assertDeb(ref.isInBounds(index),
+            "Grid index out of bounds for particle position [" << pos.x << ", " << pos.y << ", "
+                                                               << pos.z << "]");
+  return (ref.isInBounds(index)) ? index : -1;
+}
+
+static inline IndexInt indexVFace(const Vec3 &pos, const MACGrid &ref)
+{
+  const Vec3i f = toVec3i(pos), c = toVec3i(pos - 0.5);
+  const IndexInt index = c.x * ref.getStrideX() + f.y * ref.getStrideY() + c.z * ref.getStrideZ();
+  assertDeb(ref.isInBounds(index),
+            "Grid index out of bounds for particle position [" << pos.x << ", " << pos.y << ", "
+                                                               << pos.z << "]");
+  return (ref.isInBounds(index)) ? index : -1;
+}
+
+static inline IndexInt indexWFace(const Vec3 &pos, const MACGrid &ref)
+{
+  const Vec3i f = toVec3i(pos), c = toVec3i(pos - 0.5);
+  const IndexInt index = c.x * ref.getStrideX() + c.y * ref.getStrideY() + f.z * ref.getStrideZ();
+  assertDeb(ref.isInBounds(index),
+            "Grid index out of bounds for particle position [" << pos.x << ", " << pos.y << ", "
+                                                               << pos.z << "]");
+  return (ref.isInBounds(index)) ? index : -1;
+}
+
+static inline IndexInt indexOffset(
+    const IndexInt gidx, const int i, const int j, const int k, const MACGrid &ref)
+{
+  const IndexInt dX[2] = {0, ref.getStrideX()};
+  const IndexInt dY[2] = {0, ref.getStrideY()};
+  const IndexInt dZ[2] = {0, ref.getStrideZ()};
+  const IndexInt index = gidx + dX[i] + dY[j] + dZ[k];
+  assertDeb(ref.isInBounds(index),
+            "Grid index out of bounds for particle position [" << pos.x << ", " << pos.y << ", "
+                                                               << pos.z << "]");
+  return (ref.isInBounds(index)) ? index : -1;
+}
+
 struct knApicMapLinearVec3ToMACGrid : public KernelBase {
   knApicMapLinearVec3ToMACGrid(const BasicParticleSystem &p,
                                MACGrid &mg,
@@ -30,7 +78,8 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
                                const ParticleDataImpl<Vec3> &cpy,
                                const ParticleDataImpl<Vec3> &cpz,
                                const ParticleDataImpl<int> *ptype,
-                               const int exclude)
+                               const int exclude,
+                               const int boundaryWidth)
       : KernelBase(p.size()),
         p(p),
         mg(mg),
@@ -40,7 +89,8 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
         cpy(cpy),
         cpz(cpz),
         ptype(ptype),
-        exclude(exclude)
+        exclude(exclude),
+        boundaryWidth(boundaryWidth)
   {
     runMessage();
     run();
@@ -54,73 +104,91 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
                  const ParticleDataImpl<Vec3> &cpy,
                  const ParticleDataImpl<Vec3> &cpz,
                  const ParticleDataImpl<int> *ptype,
-                 const int exclude)
+                 const int exclude,
+                 const int boundaryWidth)
   {
     if (!p.isActive(idx) || (ptype && ((*ptype)[idx] & exclude)))
       return;
-    const IndexInt dX[2] = {0, vg.getStrideX()};
-    const IndexInt dY[2] = {0, vg.getStrideY()};
-    const IndexInt dZ[2] = {0, vg.getStrideZ()};
-
-    const Vec3 &pos = p[idx].pos, &vel = vp[idx];
-    const IndexInt fi = static_cast<IndexInt>(pos.x), fj = static_cast<IndexInt>(pos.y),
-                   fk = static_cast<IndexInt>(pos.z);
-    const IndexInt ci = static_cast<IndexInt>(pos.x - 0.5),
-                   cj = static_cast<IndexInt>(pos.y - 0.5),
-                   ck = static_cast<IndexInt>(pos.z - 0.5);
-    const Real wfi = clamp(pos.x - fi, Real(0), Real(1)),
-               wfj = clamp(pos.y - fj, Real(0), Real(1)),
-               wfk = clamp(pos.z - fk, Real(0), Real(1));
-    const Real wci = clamp(Real(pos.x - ci - 0.5), Real(0), Real(1)),
-               wcj = clamp(Real(pos.y - cj - 0.5), Real(0), Real(1)),
-               wck = clamp(Real(pos.z - ck - 0.5), Real(0), Real(1));
-    // TODO: check index for safety
+    if (!vg.isInBounds(p.getPos(idx), boundaryWidth)) {
+      debMsg("Skipping particle at index " << idx
+                                           << ". Is out of bounds and cannot be applied to grid.",
+             1);
+      return;
+    }
+
+    const Vec3 &pos = p.getPos(idx), &vel = vp[idx];
+    const Vec3i f = toVec3i(pos);
+    const Vec3i c = toVec3i(pos - 0.5);
+    const Vec3 wf = clamp(pos - toVec3(f), Vec3(0.), Vec3(1.));
+    const Vec3 wc = clamp(pos - toVec3(c) - 0.5, Vec3(0.), Vec3(1.));
+
     {  // u-face
-      const IndexInt gidx = fi * dX[1] + cj * dY[1] + ck * dZ[1];
-      const Vec3 gpos(fi, cj + 0.5, ck + 0.5);
-      const Real wi[2] = {Real(1) - wfi, wfi};
-      const Real wj[2] = {Real(1) - wcj, wcj};
-      const Real wk[2] = {Real(1) - wck, wck};
-      for (int i = 0; i < 2; ++i)
-        for (int j = 0; j < 2; ++j)
-          for (int k = 0; k < 2; ++k) {
-            const Real w = wi[i] * wj[j] * wk[k];
-            mg[gidx + dX[i] + dY[j] + dZ[k]].x += w;
-            vg[gidx + dX[i] + dY[j] + dZ[k]].x += w * vel.x;
-            vg[gidx + dX[i] + dY[j] + dZ[k]].x += w * dot(cpx[idx], gpos + Vec3(i, j, k) - pos);
-          }
+      const IndexInt gidx = indexUFace(pos, vg);
+      if (gidx < 0)
+        return;  // debug will fail
+
+      const Vec3 gpos(f.x, c.y + 0.5, c.z + 0.5);
+      const Real wi[2] = {Real(1) - wf.x, wf.x};
+      const Real wj[2] = {Real(1) - wc.y, wc.y};
+      const Real wk[2] = {Real(1) - wc.z, wc.z};
+
+      FOR_INT_IJK(2)
+      {
+        const Real w = wi[i] * wj[j] * wk[k];
+        const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+        if (vidx < 0)
+          continue;  // debug will fail
+
+        mg[vidx].x += w;
+        vg[vidx].x += w * vel.x;
+        vg[vidx].x += w * dot(cpx[idx], gpos + Vec3(i, j, k) - pos);
+      }
     }
     {  // v-face
-      const IndexInt gidx = ci * dX[1] + fj * dY[1] + ck * dZ[1];
-      const Vec3 gpos(ci + 0.5, fj, ck + 0.5);
-      const Real wi[2] = {Real(1) - wci, wci};
-      const Real wj[2] = {Real(1) - wfj, wfj};
-      const Real wk[2] = {Real(1) - wck, wck};
-      for (int i = 0; i < 2; ++i)
-        for (int j = 0; j < 2; ++j)
-          for (int k = 0; k < 2; ++k) {
-            const Real w = wi[i] * wj[j] * wk[k];
-            mg[gidx + dX[i] + dY[j] + dZ[k]].y += w;
-            vg[gidx + dX[i] + dY[j] + dZ[k]].y += w * vel.y;
-            vg[gidx + dX[i] + dY[j] + dZ[k]].y += w * dot(cpy[idx], gpos + Vec3(i, j, k) - pos);
-          }
+      const IndexInt gidx = indexVFace(pos, vg);
+      if (gidx < 0)
+        return;
+
+      const Vec3 gpos(c.x + 0.5, f.y, c.z + 0.5);
+      const Real wi[2] = {Real(1) - wc.x, wc.x};
+      const Real wj[2] = {Real(1) - wf.y, wf.y};
+      const Real wk[2] = {Real(1) - wc.z, wc.z};
+
+      FOR_INT_IJK(2)
+      {
+        const Real w = wi[i] * wj[j] * wk[k];
+        const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+        if (vidx < 0)
+          continue;
+
+        mg[vidx].y += w;
+        vg[vidx].y += w * vel.y;
+        vg[vidx].y += w * dot(cpy[idx], gpos + Vec3(i, j, k) - pos);
+      }
     }
     if (!vg.is3D())
       return;
     {  // w-face
-      const IndexInt gidx = ci * dX[1] + cj * dY[1] + fk * dZ[1];
-      const Vec3 gpos(ci + 0.5, cj + 0.5, fk);
-      const Real wi[2] = {Real(1) - wci, wci};
-      const Real wj[2] = {Real(1) - wcj, wcj};
-      const Real wk[2] = {Real(1) - wfk, wfk};
-      for (int i = 0; i < 2; ++i)
-        for (int j = 0; j < 2; ++j)
-          for (int k = 0; k < 2; ++k) {
-            const Real w = wi[i] * wj[j] * wk[k];
-            mg[gidx + dX[i] + dY[j] + dZ[k]].z += w;
-            vg[gidx + dX[i] + dY[j] + dZ[k]].z += w * vel.z;
-            vg[gidx + dX[i] + dY[j] + dZ[k]].z += w * dot(cpz[idx], gpos + Vec3(i, j, k) - pos);
-          }
+      const IndexInt gidx = indexWFace(pos, vg);
+      if (gidx < 0)
+        return;
+
+      const Vec3 gpos(c.x + 0.5, c.y + 0.5, f.z);
+      const Real wi[2] = {Real(1) - wc.x, wc.x};
+      const Real wj[2] = {Real(1) - wc.y, wc.y};
+      const Real wk[2] = {Real(1) - wf.z, wf.z};
+
+      FOR_INT_IJK(2)
+      {
+        const Real w = wi[i] * wj[j] * wk[k];
+        const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+        if (vidx < 0)
+          continue;
+
+        mg[vidx].z += w;
+        vg[vidx].z += w * vel.z;
+        vg[vidx].z += w * dot(cpz[idx], gpos + Vec3(i, j, k) - pos);
+      }
     }
   }
   inline const BasicParticleSystem &getArg0()
@@ -168,6 +236,11 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
     return exclude;
   }
   typedef int type8;
+  inline const int &getArg9()
+  {
+    return boundaryWidth;
+  }
+  typedef int type9;
   void runMessage()
   {
     debMsg("Executing kernel knApicMapLinearVec3ToMACGrid ", 3);
@@ -179,7 +252,7 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
   {
     const IndexInt _sz = size;
     for (IndexInt i = 0; i < _sz; i++)
-      op(i, p, mg, vg, vp, cpx, cpy, cpz, ptype, exclude);
+      op(i, p, mg, vg, vp, cpx, cpy, cpz, ptype, exclude, boundaryWidth);
   }
   const BasicParticleSystem &p;
   MACGrid &mg;
@@ -190,6 +263,7 @@ struct knApicMapLinearVec3ToMACGrid : public KernelBase {
   const ParticleDataImpl<Vec3> &cpz;
   const ParticleDataImpl<int> *ptype;
   const int exclude;
+  const int boundaryWidth;
 };
 
 void apicMapPartsToMAC(const FlagGrid &flags,
@@ -201,23 +275,22 @@ void apicMapPartsToMAC(const FlagGrid &flags,
                        const ParticleDataImpl<Vec3> &cpz,
                        MACGrid *mass = NULL,
                        const ParticleDataImpl<int> *ptype = NULL,
-                       const int exclude = 0)
+                       const int exclude = 0,
+                       const int boundaryWidth = 0)
 {
-  // affine map
-  // let's assume that the particle mass is constant, 1.0
-  const bool freeMass = !mass;
-  if (!mass)
-    mass = new MACGrid(flags.getParent());
-  else
-    mass->clear();
+  // affine map: let's assume that the particle mass is constant, 1.0
+  if (!mass) {
+    MACGrid tmpmass(vel.getParent());
+    mass = &tmpmass;
+  }
 
+  mass->clear();
   vel.clear();
-  knApicMapLinearVec3ToMACGrid(parts, *mass, vel, partVel, cpx, cpy, cpz, ptype, exclude);
+
+  knApicMapLinearVec3ToMACGrid(
+      parts, *mass, vel, partVel, cpx, cpy, cpz, ptype, exclude, boundaryWidth);
   mass->stomp(VECTOR_EPSILON);
   vel.safeDivide(*mass);
-
-  if (freeMass)
-    delete mass;
 }
 static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
 {
@@ -241,8 +314,10 @@ static PyObject *_W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
           "ptype", 8, NULL, &_lock);
       const int exclude = _args.getOpt<int>("exclude", 9, 0, &_lock);
+      const int boundaryWidth = _args.getOpt<int>("boundaryWidth", 10, 0, &_lock);
       _retval = getPyNone();
-      apicMapPartsToMAC(flags, vel, parts, partVel, cpx, cpy, cpz, mass, ptype, exclude);
+      apicMapPartsToMAC(
+          flags, vel, parts, partVel, cpx, cpy, cpz, mass, ptype, exclude, boundaryWidth);
       _args.check();
     }
     pbFinalizePlugin(parent, "apicMapPartsToMAC", !noTiming);
@@ -270,7 +345,8 @@ struct knApicMapLinearMACGridToVec3 : public KernelBase {
                                const MACGrid &vg,
                                const FlagGrid &flags,
                                const ParticleDataImpl<int> *ptype,
-                               const int exclude)
+                               const int exclude,
+                               const int boundaryWidth)
       : KernelBase(vp.size()),
         vp(vp),
         cpx(cpx),
@@ -280,7 +356,8 @@ struct knApicMapLinearMACGridToVec3 : public KernelBase {
         vg(vg),
         flags(flags),
         ptype(ptype),
-        exclude(exclude)
+        exclude(exclude),
+        boundaryWidth(boundaryWidth)
   {
     runMessage();
     run();
@@ -294,78 +371,94 @@ struct knApicMapLinearMACGridToVec3 : public KernelBase {
                  const MACGrid &vg,
                  const FlagGrid &flags,
                  const ParticleDataImpl<int> *ptype,
-                 const int exclude) const
+                 const int exclude,
+                 const int boundaryWidth) const
   {
     if (!p.isActive(idx) || (ptype && ((*ptype)[idx] & exclude)))
       return;
+    if (!vg.isInBounds(p.getPos(idx), boundaryWidth)) {
+      debMsg("Skipping particle at index " << idx
+                                           << ". Is out of bounds and cannot get value from grid.",
+             1);
+      return;
+    }
 
     vp[idx] = cpx[idx] = cpy[idx] = cpz[idx] = Vec3(Real(0));
-    const IndexInt dX[2] = {0, vg.getStrideX()}, dY[2] = {0, vg.getStrideY()},
-                   dZ[2] = {0, vg.getStrideZ()};
     const Real gw[2] = {-Real(1), Real(1)};
 
-    const Vec3 &pos = p[idx].pos;
-    const IndexInt fi = static_cast<IndexInt>(pos.x), fj = static_cast<IndexInt>(pos.y),
-                   fk = static_cast<IndexInt>(pos.z);
-    const IndexInt ci = static_cast<IndexInt>(pos.x - 0.5),
-                   cj = static_cast<IndexInt>(pos.y - 0.5),
-                   ck = static_cast<IndexInt>(pos.z - 0.5);
-    const Real wfi = clamp(pos.x - fi, Real(0), Real(1)),
-               wfj = clamp(pos.y - fj, Real(0), Real(1)),
-               wfk = clamp(pos.z - fk, Real(0), Real(1));
-    const Real wci = clamp(Real(pos.x - ci - 0.5), Real(0), Real(1)),
-               wcj = clamp(Real(pos.y - cj - 0.5), Real(0), Real(1)),
-               wck = clamp(Real(pos.z - ck - 0.5), Real(0), Real(1));
-    // TODO: check index for safety
-    {  // u
-      const IndexInt gidx = fi * dX[1] + cj * dY[1] + ck * dZ[1];
-      const Real wx[2] = {Real(1) - wfi, wfi};
-      const Real wy[2] = {Real(1) - wcj, wcj};
-      const Real wz[2] = {Real(1) - wck, wck};
-      for (int i = 0; i < 2; ++i)
-        for (int j = 0; j < 2; ++j)
-          for (int k = 0; k < 2; ++k) {
-            const IndexInt vidx = gidx + dX[i] + dY[j] + dZ[k];
-            Real vgx = vg[vidx].x;
-            vp[idx].x += wx[i] * wy[j] * wz[k] * vgx;
-            cpx[idx].x += gw[i] * wy[j] * wz[k] * vgx;
-            cpx[idx].y += wx[i] * gw[j] * wz[k] * vgx;
-            cpx[idx].z += wx[i] * wy[j] * gw[k] * vgx;
-          }
+    const Vec3 &pos = p.getPos(idx);
+    const Vec3i f = toVec3i(pos);
+    const Vec3i c = toVec3i(pos - 0.5);
+    const Vec3 wf = clamp(pos - toVec3(f), Vec3(0.), Vec3(1.));
+    const Vec3 wc = clamp(pos - toVec3(c) - 0.5, Vec3(0.), Vec3(1.));
+
+    {  // u-face
+      const IndexInt gidx = indexUFace(pos, vg);
+      if (gidx < 0)
+        return;  // debug will fail
+
+      const Real wx[2] = {Real(1) - wf.x, wf.x};
+      const Real wy[2] = {Real(1) - wc.y, wc.y};
+      const Real wz[2] = {Real(1) - wc.z, wc.z};
+
+      FOR_INT_IJK(2)
+      {
+        const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+        if (vidx < 0)
+          continue;  // debug will fail
+
+        const Real vgx = vg[vidx].x;
+        vp[idx].x += wx[i] * wy[j] * wz[k] * vgx;
+        cpx[idx].x += gw[i] * wy[j] * wz[k] * vgx;
+        cpx[idx].y += wx[i] * gw[j] * wz[k] * vgx;
+        cpx[idx].z += wx[i] * wy[j] * gw[k] * vgx;
+      }
     }
-    {  // v
-      const IndexInt gidx = ci * dX[1] + fj * dY[1] + ck * dZ[1];
-      const Real wx[2] = {Real(1) - wci, wci};
-      const Real wy[2] = {Real(1) - wfj, wfj};
-      const Real wz[2] = {Real(1) - wck, wck};
-      for (int i = 0; i < 2; ++i)
-        for (int j = 0; j < 2; ++j)
-          for (int k = 0; k < 2; ++k) {
-            const IndexInt vidx = gidx + dX[i] + dY[j] + dZ[k];
-            Real vgy = vg[vidx].y;
-            vp[idx].y += wx[i] * wy[j] * wz[k] * vgy;
-            cpy[idx].x += gw[i] * wy[j] * wz[k] * vgy;
-            cpy[idx].y += wx[i] * gw[j] * wz[k] * vgy;
-            cpy[idx].z += wx[i] * wy[j] * gw[k] * vgy;
-          }
+    {  // v-face
+      const IndexInt gidx = indexVFace(pos, vg);
+      if (gidx < 0)
+        return;
+
+      const Real wx[2] = {Real(1) - wc.x, wc.x};
+      const Real wy[2] = {Real(1) - wf.y, wf.y};
+      const Real wz[2] = {Real(1) - wc.z, wc.z};
+
+      FOR_INT_IJK(2)
+      {
+        const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+        if (vidx < 0)
+          continue;
+
+        const Real vgy = vg[vidx].y;
+        vp[idx].y += wx[i] * wy[j] * wz[k] * vgy;
+        cpy[idx].x += gw[i] * wy[j] * wz[k] * vgy;
+        cpy[idx].y += wx[i] * gw[j] * wz[k] * vgy;
+        cpy[idx].z += wx[i] * wy[j] * gw[k] * vgy;
+      }
     }
     if (!vg.is3D())
       return;
-    {  // w
-      const IndexInt gidx = ci * dX[1] + cj * dY[1] + fk * dZ[1];
-      const Real wx[2] = {Real(1) - wci, wci};
-      const Real wy[2] = {Real(1) - wcj, wcj};
-      const Real wz[2] = {Real(1) - wfk, wfk};
-      for (int i = 0; i < 2; ++i)
-        for (int j = 0; j < 2; ++j)
-          for (int k = 0; k < 2; ++k) {
-            const IndexInt vidx = gidx + dX[i] + dY[j] + dZ[k];
-            Real vgz = vg[vidx].z;
-            vp[idx].z += wx[i] * wy[j] * wz[k] * vgz;
-            cpz[idx].x += gw[i] * wy[j] * wz[k] * vgz;
-            cpz[idx].y += wx[i] * gw[j] * wz[k] * vgz;
-            cpz[idx].z += wx[i] * wy[j] * gw[k] * vgz;
-          }
+    {  // w-face
+      const IndexInt gidx = indexWFace(pos, vg);
+      if (gidx < 0)
+        return;
+
+      const Real wx[2] = {Real(1) - wc.x, wc.x};
+      const Real wy[2] = {Real(1) - wc.y, wc.y};
+      const Real wz[2] = {Real(1) - wf.z, wf.z};
+
+      FOR_INT_IJK(2)
+      {
+        const IndexInt vidx = indexOffset(gidx, i, j, k, vg);
+        if (vidx < 0)
+          continue;
+
+        const Real vgz = vg[vidx].z;
+        vp[idx].z += wx[i] * wy[j] * wz[k] * vgz;
+        cpz[idx].x += gw[i] * wy[j] * wz[k] * vgz;
+        cpz[idx].y += wx[i] * gw[j] * wz[k] * vgz;
+        cpz[idx].z += wx[i] * wy[j] * gw[k] * vgz;
+      }
     }
   }
   inline ParticleDataImpl<Vec3> &getArg0()
@@ -413,6 +506,11 @@ struct knApicMapLinearMACGridToVec3 : public KernelBase {
     return exclude;
   }
   typedef int type8;
+  inline const int &getArg9()
+  {
+    return boundaryWidth;
+  }
+  typedef int type9;
   void runMessage()
   {
     debMsg("Executing kernel knApicMapLinearMACGridToVec3 ", 3);
@@ -423,7 +521,7 @@ struct knApicMapLinearMACGridToVec3 : public KernelBase {
   void operator()(const tbb::blocked_range<IndexInt> &__r) const
   {
     for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++)
-      op(idx, vp, cpx, cpy, cpz, p, vg, flags, ptype, exclude);
+      op(idx, vp, cpx, cpy, cpz, p, vg, flags, ptype, exclude, boundaryWidth);
   }
   void run()
   {
@@ -438,6 +536,7 @@ struct knApicMapLinearMACGridToVec3 : public KernelBase {
   const FlagGrid &flags;
   const ParticleDataImpl<int> *ptype;
   const int exclude;
+  const int boundaryWidth;
 };
 
 void apicMapMACGridToParts(ParticleDataImpl<Vec3> &partVel,
@@ -448,9 +547,11 @@ void apicMapMACGridToParts(ParticleDataImpl<Vec3> &partVel,
                            const MACGrid &vel,
                            const FlagGrid &flags,
                            const ParticleDataImpl<int> *ptype = NULL,
-                           const int exclude = 0)
+                           const int exclude = 0,
+                           const int boundaryWidth = 0)
 {
-  knApicMapLinearMACGridToVec3(partVel, cpx, cpy, cpz, parts, vel, flags, ptype, exclude);
+  knApicMapLinearMACGridToVec3(
+      partVel, cpx, cpy, cpz, parts, vel, flags, ptype, exclude, boundaryWidth);
 }
 static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
 {
@@ -473,8 +574,10 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       const ParticleDataImpl<int> *ptype = _args.getPtrOpt<ParticleDataImpl<int>>(
           "ptype", 7, NULL, &_lock);
       const int exclude = _args.getOpt<int>("exclude", 8, 0, &_lock);
+      const int boundaryWidth = _args.getOpt<int>("boundaryWidth", 9, 0, &_lock);
       _retval = getPyNone();
-      apicMapMACGridToParts(partVel, cpx, cpy, cpz, parts, vel, flags, ptype, exclude);
+      apicMapMACGridToParts(
+          partVel, cpx, cpy, cpz, parts, vel, flags, ptype, exclude, boundaryWidth);
       _args.check();
     }
     pbFinalizePlugin(parent, "apicMapMACGridToParts", !noTiming);
diff --git a/intern/mantaflow/intern/MANTA_main.cpp b/intern/mantaflow/intern/MANTA_main.cpp
index e15fe809b2d..729cb151a39 100644
--- a/intern/mantaflow/intern/MANTA_main.cpp
+++ b/intern/mantaflow/intern/MANTA_main.cpp
@@ -697,12 +697,6 @@ void MANTA::initializeRNAMap(FluidModifierData *fmd)
   if ((fds->border_collisions & FLUID_DOMAIN_BORDER_TOP) == 0)
     borderCollisions += "Z";
 
-  string simulationMethod = "";
-  if (fds->simulation_method & FLUID_DOMAIN_METHOD_FLIP)
-    simulationMethod += "'FLIP'";
-  else if (fds->simulation_method & FLUID_DOMAIN_METHOD_APIC)
-    simulationMethod += "'APIC'";
-
   string particleTypesStr = "";
   if (fds->particle_type & FLUID_DOMAIN_PARTICLE_SPRAY)
     particleTypesStr += "PtypeSpray";
@@ -837,7 +831,7 @@ void MANTA::initializeRNAMap(FluidModifierData *fmd)
   mRNAMap["CACHE_MESH_FORMAT"] = getCacheFileEnding(fds->cache_mesh_format);
   mRNAMap["CACHE_NOISE_FORMAT"] = getCacheFileEnding(fds->cache_noise_format);
   mRNAMap["CACHE_PARTICLE_FORMAT"] = getCacheFileEnding(fds->cache_particle_format);
-  mRNAMap["SIMULATION_METHOD"] = simulationMethod;
+  mRNAMap["USING_APIC"] = getBooleanString(fds->simulation_method == FLUID_DOMAIN_METHOD_APIC);
   mRNAMap["FLIP_RATIO"] = to_string(fds->flip_ratio);
   mRNAMap["PARTICLE_RANDOMNESS"] = to_string(fds->particle_randomness);
   mRNAMap["PARTICLE_NUMBER"] = to_string(fds->particle_number);
diff --git a/intern/mantaflow/intern/strings/liquid_script.h b/intern/mantaflow/intern/strings/liquid_script.h
index 7622fc5a602..8e0d113d680 100644
--- a/intern/mantaflow/intern/strings/liquid_script.h
+++ b/intern/mantaflow/intern/strings/liquid_script.h
@@ -50,7 +50,8 @@ smoothenPos_s$ID$      = $MESH_SMOOTHEN_POS$\n\
 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";
+maxSysParticles_s$ID$  = $PP_PARTICLE_MAXIMUM$\n\
+using_apic_s$ID$       = $USING_APIC$\n";
 
 const std::string liquid_variables_particles =
     "\n\
@@ -91,6 +92,14 @@ curvature_s$ID$  = None\n\
 pp_s$ID$         = s$ID$.create(BasicParticleSystem, name='$NAME_PARTS$')\n\
 pVel_pp$ID$      = pp_s$ID$.create(PdataVec3, name='$NAME_PARTSVELOCITY$')\n\
 \n\
+pCx_pp$ID$       = None\n\
+pCy_pp$ID$       = None\n\
+pCz_pp$ID$       = None\n\
+if using_apic_s$ID$:\n\
+    pCx_pp$ID$   = pp_s$ID$.create(PdataVec3)\n\
+    pCy_pp$ID$   = pp_s$ID$.create(PdataVec3)\n\
+    pCz_pp$ID$   = pp_s$ID$.create(PdataVec3)\n\
+\n\
 # Acceleration data for particle nbs\n\
 pindex_s$ID$     = s$ID$.create(ParticleIndexSystem, name='$NAME_PINDEX$')\n\
 gpi_s$ID$        = s$ID$.create(IntGrid, name='$NAME_GPI$')\n\
@@ -275,8 +284,12 @@ def liquid_step_$ID$():\n\
         resetOutflow(flags=flags_s$ID$, phi=phi_s$ID$, parts=pp_s$ID$, index=gpi_s$ID$, indexSys=pindex_s$ID$)\n\
     flags_s$ID$.updateFromLevelset(phi_s$ID$)\n\
     \n\
-    # combine particles velocities with advected grid velocities\n\
-    mapPartsToMAC(vel=velParts_s$ID$, flags=flags_s$ID$, velOld=velOld_s$ID$, parts=pp_s$ID$, partVel=pVel_pp$ID$, weight=mapWeights_s$ID$)\n\
+    # combine particle velocities with advected grid velocities\n\
+    if using_apic_s$ID$:\n\
+        apicMapPartsToMAC(flags=flags_s$ID$, vel=vel_s$ID$, parts=pp_s$ID$, partVel=pVel_pp$ID$, cpx=pCx_pp$ID$, cpy=pCy_pp$ID$, cpz=pCz_pp$ID$)\n\
+    else:\n\
+        mapPartsToMAC(vel=velParts_s$ID$, flags=flags_s$ID$, velOld=velOld_s$ID$, parts=pp_s$ID$, partVel=pVel_pp$ID$, weight=mapWeights_s$ID$)\n\
+    \n\
     extrapolateMACFromWeight(vel=velParts_s$ID$, distance=2, weight=mapWeights_s$ID$)\n\
     combineGridVel(vel=velParts_s$ID$, weight=mapWeights_s$ID$, combineVel=vel_s$ID$, phi=phi_s$ID$, narrowBand=combineBandWidth_s$ID$, thresh=0)\n\
     velOld_s$ID$.copyFrom(vel_s$ID$)\n\
@@ -319,7 +332,11 @@ def liquid_step_$ID$():\n\
     # set source grids for resampling, used in adjustNumber!\n\
     pVel_pp$ID$.setSource(grid=vel_s$ID$, isMAC=True)\n\
     adjustNumber(parts=pp_s$ID$, vel=vel_s$ID$, flags=flags_s$ID$, minParticles=minParticles_s$ID$, maxParticles=maxParticles_s$ID$, phi=phi_s$ID$, exclude=phiObs_s$ID$, radiusFactor=radiusFactor_s$ID$, narrowBand=adjustedNarrowBandWidth_s$ID$)\n\
-    flipVelocityUpdate(vel=vel_s$ID$, velOld=velOld_s$ID$, flags=flags_s$ID$, parts=pp_s$ID$, partVel=pVel_pp$ID$, flipRatio=flipRatio_s$ID$)\n";
+    \n\
+    if using_apic_s$ID$:\n\
+        apicMapMACGridToParts(partVel=pVel_pp$ID$, cpx=pCx_pp$ID$, cpy=pCy_pp$ID$, cpz=pCz_pp$ID$, parts=pp_s$ID$, vel=vel_s$ID$, flags=flags_s$ID$)\n\
+    else:\n\
+        flipVelocityUpdate(vel=vel_s$ID$, velOld=velOld_s$ID$, flags=flags_s$ID$, parts=pp_s$ID$, partVel=pVel_pp$ID$, flipRatio=flipRatio_s$ID$)\n";
 
 const std::string liquid_step_mesh =
     "\n\
diff --git a/release/scripts/startup/bl_ui/properties_physics_fluid.py b/release/scripts/startup/bl_ui/properties_physics_fluid.py
index 15c5e70d8b2..3afe7a47028 100644
--- a/release/scripts/startup/bl_ui/properties_physics_fluid.py
+++ b/release/scripts/startup/bl_ui/properties_physics_fluid.py
@@ -486,7 +486,8 @@ class PHYSICS_PT_liquid(PhysicButtonsPanel, Panel):
 
         col = flow.column()
         col.prop(domain, "simulation_method", expand=False)
-        col.prop(domain, "flip_ratio", text="FLIP Ratio")
+        if domain.simulation_method == 'FLIP':
+            col.prop(domain, "flip_ratio", text="FLIP Ratio")
         col.prop(domain, "sys_particle_maximum", text="System Maximum")
         col = col.column(align=True)
         col.prop(domain, "particle_radius", text="Particle Radius")
diff --git a/source/blender/blenkernel/intern/fluid.c b/source/blender/blenkernel/intern/fluid.c
index fa1b1997625..af4b98e1c8d 100644
--- a/source/blender/blenkernel/intern/fluid.c
+++ b/source/blender/blenkernel/intern/fluid.c
@@ -4993,6 +4993,7 @@ void BKE_fluid_modifier_copy(const struct FluidModifierData *fmd,
     tfds->fractions_threshold = fds->fractions_threshold;
     tfds->fractions_distance = fds->fractions_distance;
     tfds->sys_particle_maximum = fds->sys_particle_maximum;
+    tfds->simulation_method = fds->simulation_method;
 
     /* diffusion options*/
     tfds->surface_tension = fds->surface_tension;
diff --git a/source/blender/makesrna/intern/rna_fluid.c b/source/blender/makesrna/intern/rna_fluid.c
index 60eed60ace1..02ae71e3b3d 100644
--- a/source/blender/makesrna/intern/rna_fluid.c
+++ b/source/blender/makesrna/intern/rna_fluid.c
@@ -1456,8 +1456,16 @@ static void rna_def_fluid_domain_settings(BlenderRNA *brna)
       {0, NULL, 0, NULL, NULL}};
 
   static EnumPropertyItem simulation_methods[] = {
-      {FLUID_DOMAIN_METHOD_FLIP, "FLIP", 0, "FLIP", "Use FLIP as the simulation method"},
-      /*{FLUID_DOMAIN_METHOD_APIC, "APIC", 0, "APIC", "Use APIC as the simulation method"},*/
+      {FLUID_DOMAIN_METHOD_FLIP,
+       "FLIP",
+       0,
+       "FLIP",
+       "Use FLIP as the simulation method (more splashy behavior)"},
+      {FLUID_DOMAIN_METHOD_APIC,
+       "APIC",
+       0,
+       "APIC",
+       "Use APIC as the simulation method (more energetic and stable behavior)"},
       {0, NULL, 0, NULL, NULL},
   };
 
@@ -1820,6 +1828,7 @@ static void rna_def_fluid_domain_settings(BlenderRNA *brna)
   RNA_def_property_enum_sdna(prop, NULL, "simulation_method");
   RNA_def_property_enum_items(prop, simulation_methods);
   RNA_def_property_ui_text(prop, "Simulation Method", "Change the underlying simulation method");
+  RNA_def_property_clear_flag(prop, PROP_ANIMATABLE);
   RNA_def_property_update(prop, NC_OBJECT | ND_MODIFIER, "rna_Fluid_domain_data_reset");
 
   prop = RNA_def_property(srna, "flip_ratio", PROP_FLOAT, PROP_NONE);