Fluid: Updated Mantaflow source files

This update introduces two improvements from the Mantaflow repository:

(1) Improved particle sampling:
- Liquid and secondary particles are sampled more predictably. With all parameters being equal, baked particles will be computed at the exact same position during every bake.
- Before, this was not guaranteed.

(2) Sparse grid caching:
- While saving grid data to disk, grids will from now on be saved in a sparse structure whenever possible (e.g. density, flame but not levelsets).
- With the sparse optimization grid cells with a value under the 'Empty Space' value (already present in domain settings) will not be cached.
- The main benefits of this optimization are: Smaller cache sizes and faster playback of simulation data in the viewport.
- This optimization works 'out-of-the-box'. There is no option in the UI to enable it.
- For now, only smoke simulation grids will take advantage of this optimization.

10 files changed, 239 insertions, 97 deletions

diff --git a/extern/mantaflow/preprocessed/fileio/iovdb.cpp b/extern/mantaflow/preprocessed/fileio/iovdb.cpp
index 287c78f0608..8aca19ade04 100644
--- a/extern/mantaflow/preprocessed/fileio/iovdb.cpp
+++ b/extern/mantaflow/preprocessed/fileio/iovdb.cpp
@@ -31,6 +31,8 @@
 #  include "openvdb/openvdb.h"
 #  include <openvdb/points/PointConversion.h>
 #  include <openvdb/points/PointCount.h>
+#  include <openvdb/tools/Clip.h>
+#  include <openvdb/tools/Dense.h>
 #endif
 
 #define POSITION_NAME "P"
@@ -45,15 +47,29 @@ namespace Manta {
 template<class GridType, class T> void importVDB(typename GridType::Ptr from, Grid<T> *to)
 {
   using ValueT = typename GridType::ValueType;
-  typename GridType::Accessor accessor = from->getAccessor();
 
-  FOR_IJK(*to)
-  {
-    openvdb::Coord xyz(i, j, k);
-    ValueT vdbValue = accessor.getValue(xyz);
-    T toMantaValue;
-    convertFrom(vdbValue, &toMantaValue);
-    to->set(i, j, k, toMantaValue);
+  // Check if current grid is to be read as a sparse grid, active voxels (only) will be copied
+  if (to->saveSparse()) {
+    to->clear();  // Ensure that destination grid is empty before writing
+    for (typename GridType::ValueOnCIter iter = from->cbeginValueOn(); iter.test(); ++iter) {
+      ValueT vdbValue = *iter;
+      openvdb::Coord coord = iter.getCoord();
+      T toMantaValue;
+      convertFrom(vdbValue, &toMantaValue);
+      to->set(coord.x(), coord.y(), coord.z(), toMantaValue);
+    }
+  }
+  // When importing all grid cells, using a grid accessor is usually faster than a value iterator
+  else {
+    typename GridType::Accessor accessor = from->getAccessor();
+    FOR_IJK(*to)
+    {
+      openvdb::Coord xyz(i, j, k);
+      ValueT vdbValue = accessor.getValue(xyz);
+      T toMantaValue;
+      convertFrom(vdbValue, &toMantaValue);
+      to->set(i, j, k, toMantaValue);
+    }
   }
 }
 
@@ -173,19 +189,43 @@ static void setGridOptions(typename GridType::Ptr grid,
   grid->setSaveFloatAsHalf(precision == PRECISION_MINI || precision == PRECISION_HALF);
 }
 
-template<class T, class GridType> typename GridType::Ptr exportVDB(Grid<T> *from)
+template<class T, class GridType>
+typename GridType::Ptr exportVDB(Grid<T> *from, float clip, openvdb::FloatGrid::Ptr clipGrid)
 {
   using ValueT = typename GridType::ValueType;
-  typename GridType::Ptr to = GridType::create();
-  typename GridType::Accessor accessor = to->getAccessor();
-
-  FOR_IJK(*from)
-  {
-    openvdb::Coord xyz(i, j, k);
-    T fromMantaValue = (*from)(i, j, k);
-    ValueT vdbValue;
-    convertTo(&vdbValue, fromMantaValue);
-    accessor.setValue(xyz, vdbValue);
+  typename GridType::Ptr to = GridType::create(ValueT(0));
+
+  // Copy data from grid by creating a vdb dense structure and then copying that into a vdb grid
+  // This is the fastest way to copy data for both dense and sparse grids -> if (true)
+  if (true) {
+    ValueT *data = (ValueT *)from->getData();
+    openvdb::math::CoordBBox bbox(
+        openvdb::Coord(0),
+        openvdb::Coord(from->getSizeX() - 1, from->getSizeY() - 1, from->getSizeZ() - 1));
+    openvdb::tools::Dense<ValueT, openvdb::tools::MemoryLayout::LayoutXYZ> dense(bbox, data);
+
+    // Trick: Set clip value to very small / negative value in order to copy all values of dense
+    // grids
+    float tmpClip = (from->saveSparse()) ? clip : -std::numeric_limits<Real>::max();
+    // Copy from dense to sparse grid structure considering clip value
+    openvdb::tools::copyFromDense(dense, *to, ValueT(tmpClip));
+
+    // If present, use clip grid to trim down current vdb grid even more
+    if (from->saveSparse() && clipGrid && !clipGrid->empty()) {
+      to = openvdb::tools::clip(*to, *clipGrid);
+    }
+  }
+  // Alternatively, reading all grid cells with an accessor (slightly slower) is possible like this
+  else {
+    typename GridType::Accessor accessor = to->getAccessor();
+    FOR_IJK(*from)
+    {
+      openvdb::Coord xyz(i, j, k);
+      T fromMantaValue = (*from)(i, j, k);
+      ValueT vdbValue;
+      convertTo(&vdbValue, fromMantaValue);
+      accessor.setValue(xyz, vdbValue);
+    }
   }
   return to;
 }
@@ -346,7 +386,9 @@ int writeObjectsVDB(const string &filename,
                     float worldSize,
                     bool skipDeletedParts,
                     int compression,
-                    int precision)
+                    int precision,
+                    float clip,
+                    const Grid<Real> *clipGrid)
 {
   openvdb::initialize();
   openvdb::io::File file(filename);
@@ -357,6 +399,19 @@ int writeObjectsVDB(const string &filename,
 
   std::vector<ParticleDataBase *> pdbBuffer;
 
+  // Convert given clip grid to vdb clip grid
+  openvdb::FloatGrid::Ptr vdbClipGrid = nullptr;
+  if (clipGrid) {
+    vdbClipGrid = openvdb::FloatGrid::create();
+    Real *data = (Real *)clipGrid->getData();
+    openvdb::math::CoordBBox bbox(openvdb::Coord(0),
+                                  openvdb::Coord(clipGrid->getSizeX() - 1,
+                                                 clipGrid->getSizeY() - 1,
+                                                 clipGrid->getSizeZ() - 1));
+    openvdb::tools::Dense<float, openvdb::tools::MemoryLayout::LayoutXYZ> dense(bbox, data);
+    openvdb::tools::copyFromDense(dense, *vdbClipGrid, clip);
+  }
+
   for (std::vector<PbClass *>::iterator iter = objects->begin(); iter != objects->end(); ++iter) {
     openvdb::GridClass gClass = openvdb::GRID_UNKNOWN;
     openvdb::GridBase::Ptr vdbGrid;
@@ -368,10 +423,14 @@ int writeObjectsVDB(const string &filename,
 
     if (GridBase *mantaGrid = dynamic_cast<GridBase *>(*iter)) {
 
+      if (clipGrid) {
+        assertMsg(clipGrid->getSize() == mantaGrid->getSize(),
+                  "writeObjectsVDB: Clip grid and exported grid must have the same size");
+      }
       if (mantaGrid->getType() & GridBase::TypeInt) {
         debMsg("Writing int grid '" << mantaGrid->getName() << "' to vdb file " << filename, 1);
         Grid<int> *mantaIntGrid = (Grid<int> *)mantaGrid;
-        vdbGrid = exportVDB<int, openvdb::Int32Grid>(mantaIntGrid);
+        vdbGrid = exportVDB<int, openvdb::Int32Grid>(mantaIntGrid, clip, vdbClipGrid);
         gridsVDB.push_back(vdbGrid);
       }
       else if (mantaGrid->getType() & GridBase::TypeReal) {
@@ -379,7 +438,9 @@ int writeObjectsVDB(const string &filename,
         gClass = (mantaGrid->getType() & GridBase::TypeLevelset) ? openvdb::GRID_LEVEL_SET :
                                                                    openvdb::GRID_FOG_VOLUME;
         Grid<Real> *mantaRealGrid = (Grid<Real> *)mantaGrid;
-        vdbGrid = exportVDB<Real, openvdb::FloatGrid>(mantaRealGrid);
+        // Only supply clip grid if real grid is not equal to the clip grid
+        openvdb::FloatGrid::Ptr tmpClipGrid = (mantaRealGrid == clipGrid) ? nullptr : vdbClipGrid;
+        vdbGrid = exportVDB<Real, openvdb::FloatGrid>(mantaRealGrid, clip, tmpClipGrid);
         gridsVDB.push_back(vdbGrid);
       }
       else if (mantaGrid->getType() & GridBase::TypeVec3) {
@@ -387,7 +448,7 @@ int writeObjectsVDB(const string &filename,
         gClass = (mantaGrid->getType() & GridBase::TypeMAC) ? openvdb::GRID_STAGGERED :
                                                               openvdb::GRID_UNKNOWN;
         Grid<Vec3> *mantaVec3Grid = (Grid<Vec3> *)mantaGrid;
-        vdbGrid = exportVDB<Vec3, openvdb::Vec3SGrid>(mantaVec3Grid);
+        vdbGrid = exportVDB<Vec3, openvdb::Vec3SGrid>(mantaVec3Grid, clip, vdbClipGrid);
         gridsVDB.push_back(vdbGrid);
       }
       else {
@@ -620,9 +681,12 @@ template void importVDB<openvdb::Int32Grid, int>(openvdb::Int32Grid::Ptr from, G
 template void importVDB<openvdb::FloatGrid, Real>(openvdb::FloatGrid::Ptr from, Grid<Real> *to);
 template void importVDB<openvdb::Vec3SGrid, Vec3>(openvdb::Vec3SGrid::Ptr from, Grid<Vec3> *to);
 
-template openvdb::Int32Grid::Ptr exportVDB<int, openvdb::Int32Grid>(Grid<int> *from);
-template openvdb::FloatGrid::Ptr exportVDB<Real, openvdb::FloatGrid>(Grid<Real> *from);
-template openvdb::Vec3SGrid::Ptr exportVDB<Vec3, openvdb::Vec3SGrid>(Grid<Vec3> *from);
+template openvdb::Int32Grid::Ptr exportVDB<int, openvdb::Int32Grid>(
+    Grid<int> *from, float clip = 1e-4, openvdb::FloatGrid::Ptr clipGrid = nullptr);
+template openvdb::FloatGrid::Ptr exportVDB<Real, openvdb::FloatGrid>(
+    Grid<Real> *from, float clip = 1e-4, openvdb::FloatGrid::Ptr clipGrid = nullptr);
+template openvdb::Vec3SGrid::Ptr exportVDB<Vec3, openvdb::Vec3SGrid>(
+    Grid<Vec3> *from, float clip = 1e-4, openvdb::FloatGrid::Ptr clipGrid = nullptr);
 
 openvdb::points::PointDataGrid::Ptr exportVDB(BasicParticleSystem *from,
                                               std::vector<ParticleDataBase *> &fromPData,
@@ -652,7 +716,9 @@ int writeObjectsVDB(const string &filename,
                     float worldSize,
                     bool skipDeletedParts,
                     int compression,
-                    int precision)
+                    int precision,
+                    float clip,
+                    const Grid<Real> *clipGrid)
 {
   errMsg("Cannot save to .vdb file. Mantaflow has not been built with OpenVDB support.");
   return 0;
diff --git a/extern/mantaflow/preprocessed/fileio/mantaio.cpp b/extern/mantaflow/preprocessed/fileio/mantaio.cpp
index 3048572261a..fe29890ec11 100644
--- a/extern/mantaflow/preprocessed/fileio/mantaio.cpp
+++ b/extern/mantaflow/preprocessed/fileio/mantaio.cpp
@@ -83,7 +83,9 @@ int save(const string &name,
          bool skipDeletedParts = false,
          int compression = COMPRESSION_ZIP,
          bool precisionHalf = true,
-         int precision = PRECISION_HALF)
+         int precision = PRECISION_HALF,
+         float clip = 1e-4,
+         const Grid<Real> *clipGrid = nullptr)
 {
 
   if (!precisionHalf) {
@@ -102,7 +104,8 @@ int save(const string &name,
   else if (ext == ".vol")
     return writeGridsVol(name, &objects);
   if (ext == ".vdb")
-    return writeObjectsVDB(name, &objects, worldSize, skipDeletedParts, compression, precision);
+    return writeObjectsVDB(
+        name, &objects, worldSize, skipDeletedParts, compression, precision, clip, clipGrid);
   else if (ext == ".npz")
     return writeGridsNumpy(name, &objects);
   else if (ext == ".txt")
@@ -129,8 +132,17 @@ static PyObject *_W_1(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
       int compression = _args.getOpt<int>("compression", 4, COMPRESSION_ZIP, &_lock);
       bool precisionHalf = _args.getOpt<bool>("precisionHalf", 5, true, &_lock);
       int precision = _args.getOpt<int>("precision", 6, PRECISION_HALF, &_lock);
-      _retval = toPy(
-          save(name, objects, worldSize, skipDeletedParts, compression, precisionHalf, precision));
+      float clip = _args.getOpt<float>("clip", 7, 1e-4, &_lock);
+      const Grid<Real> *clipGrid = _args.getPtrOpt<Grid<Real>>("clipGrid", 8, nullptr, &_lock);
+      _retval = toPy(save(name,
+                          objects,
+                          worldSize,
+                          skipDeletedParts,
+                          compression,
+                          precisionHalf,
+                          precision,
+                          clip,
+                          clipGrid));
       _args.check();
     }
     pbFinalizePlugin(parent, "save", !noTiming);
diff --git a/extern/mantaflow/preprocessed/fileio/mantaio.h b/extern/mantaflow/preprocessed/fileio/mantaio.h
index d49f4b1369d..088d43556e1 100644
--- a/extern/mantaflow/preprocessed/fileio/mantaio.h
+++ b/extern/mantaflow/preprocessed/fileio/mantaio.h
@@ -75,7 +75,9 @@ int writeObjectsVDB(const std::string &filename,
                     float scale = 1.0,
                     bool skipDeletedParts = false,
                     int compression = COMPRESSION_ZIP,
-                    int precision = PRECISION_HALF);
+                    int precision = PRECISION_HALF,
+                    float clip = 1e-4,
+                    const Grid<Real> *clipGrid = nullptr);
 int readObjectsVDB(const std::string &filename,
                    std::vector<PbClass *> *objects,
                    float scale = 1.0);
diff --git a/extern/mantaflow/preprocessed/gitinfo.h b/extern/mantaflow/preprocessed/gitinfo.h
index c7e62bf1248..d7a69564a3b 100644
--- a/extern/mantaflow/preprocessed/gitinfo.h
+++ b/extern/mantaflow/preprocessed/gitinfo.h
@@ -1,3 +1,3 @@
 
 
-#define MANTA_GIT_VERSION "commit 73990a8a5b876e2b136a646258f714c08c5828da"
+#define MANTA_GIT_VERSION "commit bb7cde47b6e04fa62815c70775dc70f02065599f"
diff --git a/extern/mantaflow/preprocessed/grid.cpp b/extern/mantaflow/preprocessed/grid.cpp
index cf8e4635462..61672129f37 100644
--- a/extern/mantaflow/preprocessed/grid.cpp
+++ b/extern/mantaflow/preprocessed/grid.cpp
@@ -60,7 +60,7 @@ template<> inline GridBase::GridType typeList<Vec3>()
 }
 
 template<class T>
-Grid<T>::Grid(FluidSolver *parent, bool show) : GridBase(parent), externalData(false)
+Grid<T>::Grid(FluidSolver *parent, bool show, bool sparse) : GridBase(parent), mExternalData(false)
 {
   mType = typeList<T>();
   mSize = parent->getGridSize();
@@ -70,22 +70,23 @@ Grid<T>::Grid(FluidSolver *parent, bool show) : GridBase(parent), externalData(f
   mDx = 1.0 / mSize.max();
   clear();
   setHidden(!show);
+
+#if OPENVDB == 1
+  mSaveSparse = sparse;
+#else
+  if (sparse)
+    debMsg("Cannot enable sparse save option without OpenVDB", 1);
+  mSaveSparse = false;
+#endif
 }
 
 template<class T>
-Grid<T>::Grid(FluidSolver *parent, T *data, bool show)
-    : GridBase(parent), mData(data), externalData(true)
+Grid<T>::Grid(FluidSolver *parent, T *data, bool show) : Grid<T>::Grid(parent, show)
 {
-  mType = typeList<T>();
-  mSize = parent->getGridSize();
-
-  mStrideZ = parent->is2D() ? 0 : (mSize.x * mSize.y);
-  mDx = 1.0 / mSize.max();
-
-  setHidden(!show);
+  mData = data;
 }
 
-template<class T> Grid<T>::Grid(const Grid<T> &a) : GridBase(a.getParent()), externalData(false)
+template<class T> Grid<T>::Grid(const Grid<T> &a) : GridBase(a.getParent()), mExternalData(false)
 {
   mSize = a.mSize;
   mType = a.mType;
@@ -98,7 +99,7 @@ template<class T> Grid<T>::Grid(const Grid<T> &a) : GridBase(a.getParent()), ext
 
 template<class T> Grid<T>::~Grid()
 {
-  if (!externalData) {
+  if (!mExternalData) {
     mParent->freeGridPointer<T>(mData);
   }
 }
@@ -114,8 +115,8 @@ template<class T> void Grid<T>::swap(Grid<T> &other)
       other.getSizeZ() != getSizeZ())
     errMsg("Grid::swap(): Grid dimensions mismatch.");
 
-  if (externalData || other.externalData)
-    errMsg("Grid::swap(): Cannot swap if one grid stores externalData.");
+  if (mExternalData || other.mExternalData)
+    errMsg("Grid::swap(): Cannot swap if one grid stores mExternalData.");
 
   T *dswap = other.mData;
   other.mData = mData;
diff --git a/extern/mantaflow/preprocessed/grid.h b/extern/mantaflow/preprocessed/grid.h
index 1f3dc2789ac..cf942a19e9a 100644
--- a/extern/mantaflow/preprocessed/grid.h
+++ b/extern/mantaflow/preprocessed/grid.h
@@ -402,7 +402,7 @@ class GridBase : public PbClass {
 template<class T> class Grid : public GridBase {
  public:
   //! init new grid, values are set to zero
-  Grid(FluidSolver *parent, bool show = true);
+  Grid(FluidSolver *parent, bool show = true, bool sparse = false);
   static int _W_10(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
   {
     PbClass *obj = Pb::objFromPy(_self);
@@ -416,7 +416,8 @@ template<class T> class Grid : public GridBase {
         ArgLocker _lock;
         FluidSolver *parent = _args.getPtr<FluidSolver>("parent", 0, &_lock);
         bool show = _args.getOpt<bool>("show", 1, true, &_lock);
-        obj = new Grid(parent, show);
+        bool sparse = _args.getOpt<bool>("sparse", 2, false, &_lock);
+        obj = new Grid(parent, show, sparse);
         obj->registerObject(_self, &_args);
         _args.check();
       }
@@ -581,6 +582,16 @@ template<class T> class Grid : public GridBase {
     DEBUG_ONLY(checkIndex(idx));
     return mData[idx];
   }
+  //! raw data access
+  inline T *getData() const
+  {
+    return mData;
+  }
+  //! query if this grid should be saved as a sparse grid
+  inline bool saveSparse()
+  {
+    return mSaveSparse;
+  }
 
   //! set data
   inline void set(int i, int j, int k, T &val)
@@ -1290,7 +1301,8 @@ template<class T> class Grid : public GridBase {
 
  protected:
   T *mData;
-  bool externalData;  // True if mData is managed outside of the Fluidsolver
+  bool mExternalData;  // True if mData is managed outside of the Fluidsolver
+  bool mSaveSparse;    // True if this grid may be cached in a sparse structure
  public:
   PbArgs _args;
 }
@@ -1302,7 +1314,8 @@ template<class T> class Grid : public GridBase {
 //! Special function for staggered grids
 class MACGrid : public Grid<Vec3> {
  public:
-  MACGrid(FluidSolver *parent, bool show = true) : Grid<Vec3>(parent, show)
+  MACGrid(FluidSolver *parent, bool show = true, bool sparse = false)
+      : Grid<Vec3>(parent, show, sparse)
   {
     mType = (GridType)(TypeMAC | TypeVec3);
   }
@@ -1319,7 +1332,8 @@ class MACGrid : public Grid<Vec3> {
         ArgLocker _lock;
         FluidSolver *parent = _args.getPtr<FluidSolver>("parent", 0, &_lock);
         bool show = _args.getOpt<bool>("show", 1, true, &_lock);
-        obj = new MACGrid(parent, show);
+        bool sparse = _args.getOpt<bool>("sparse", 2, false, &_lock);
+        obj = new MACGrid(parent, show, sparse);
         obj->registerObject(_self, &_args);
         _args.check();
       }
@@ -1425,7 +1439,8 @@ class MACGrid : public Grid<Vec3> {
 //! Special functions for FlagGrid
 class FlagGrid : public Grid<int> {
  public:
-  FlagGrid(FluidSolver *parent, int dim = 3, bool show = true) : Grid<int>(parent, show)
+  FlagGrid(FluidSolver *parent, int dim = 3, bool show = true, bool sparse = false)
+      : Grid<int>(parent, show, sparse)
   {
     mType = (GridType)(TypeFlags | TypeInt);
   }
@@ -1443,7 +1458,8 @@ class FlagGrid : public Grid<int> {
         FluidSolver *parent = _args.getPtr<FluidSolver>("parent", 0, &_lock);
         int dim = _args.getOpt<int>("dim", 1, 3, &_lock);
         bool show = _args.getOpt<bool>("show", 2, true, &_lock);
-        obj = new FlagGrid(parent, dim, show);
+        bool sparse = _args.getOpt<bool>("sparse", 3, false, &_lock);
+        obj = new FlagGrid(parent, dim, show, sparse);
         obj->registerObject(_self, &_args);
         _args.check();
       }
diff --git a/extern/mantaflow/preprocessed/particle.cpp b/extern/mantaflow/preprocessed/particle.cpp
index 05c561e2c60..ad1c344d307 100644
--- a/extern/mantaflow/preprocessed/particle.cpp
+++ b/extern/mantaflow/preprocessed/particle.cpp
@@ -28,9 +28,15 @@
 using namespace std;
 namespace Manta {
 
-ParticleBase::ParticleBase(FluidSolver *parent)
-    : PbClass(parent), mMaxParticles(0), mAllowCompress(true), mFreePdata(false)
+int ParticleBase::globalSeed = 9832;
+
+ParticleBase::ParticleBase(FluidSolver *parent, int fixedSeed)
+    : PbClass(parent), mMaxParticles(0), mAllowCompress(true), mFreePdata(false), mSeed(fixedSeed)
 {
+  // use global random seed if none is given
+  if (fixedSeed == -1) {
+    mSeed = globalSeed;
+  }
 }
 
 ParticleBase::~ParticleBase()
diff --git a/extern/mantaflow/preprocessed/particle.h b/extern/mantaflow/preprocessed/particle.h
index 7fed33c737d..7fcc7e5ca32 100644
--- a/extern/mantaflow/preprocessed/particle.h
+++ b/extern/mantaflow/preprocessed/particle.h
@@ -47,7 +47,7 @@ class ParticleBase : public PbClass {
     PINVALID = (1 << 30),  // unused
   };
 
-  ParticleBase(FluidSolver *parent);
+  ParticleBase(FluidSolver *parent, int fixedSeed = -1);
   static int _W_0(PyObject *_self, PyObject *_linargs, PyObject *_kwds)
   {
     PbClass *obj = Pb::objFromPy(_self);
@@ -60,7 +60,8 @@ class ParticleBase : public PbClass {
       {
         ArgLocker _lock;
         FluidSolver *parent = _args.getPtr<FluidSolver>("parent", 0, &_lock);
-        obj = new ParticleBase(parent);
+        int fixedSeed = _args.getOpt<int>("fixedSeed", 1, -1, &_lock);
+        obj = new ParticleBase(parent, fixedSeed);
         obj->registerObject(_self, &_args);
         _args.check();
       }
@@ -111,6 +112,11 @@ class ParticleBase : public PbClass {
     return;
   }
 
+  inline int getSeed()
+  {
+    return mSeed;
+  }
+
   //! particle data functions
 
   //! create a particle data object
@@ -192,9 +198,13 @@ class ParticleBase : public PbClass {
   //! per particle)
   std::vector<ParticleDataImpl<Real> *> mPdataReal;
   std::vector<ParticleDataImpl<Vec3> *> mPdataVec3;
-  std::vector<ParticleDataImpl<int> *>
-      mPdataInt;  //! indicate that pdata of this particle system is copied, and needs to be freed
+  std::vector<ParticleDataImpl<int> *> mPdataInt;
+  //! indicate that pdata of this particle system is copied, and needs to be freed
   bool mFreePdata;
+
+  //! custom seed for particle systems, used by plugins
+  int mSeed;  //! fix global random seed storage, used mainly by functions in this class
+  static int globalSeed;
  public:
   PbArgs _args;
 }
@@ -2285,15 +2295,14 @@ void ParticleSystem<S>::advectInGrid(const FlagGrid &flags,
 }
 
 template<class S> struct KnProjectParticles : public KernelBase {
-  KnProjectParticles(ParticleSystem<S> &part, Grid<Vec3> &gradient)
-      : KernelBase(part.size()), part(part), gradient(gradient)
+  KnProjectParticles(ParticleSystem<S> &part, Grid<Vec3> &gradient, RandomStream &rand)
+      : KernelBase(part.size()), part(part), gradient(gradient), rand(rand)
   {
     runMessage();
     run();
   }
-  inline void op(IndexInt idx, ParticleSystem<S> &part, Grid<Vec3> &gradient)
+  inline void op(IndexInt idx, ParticleSystem<S> &part, Grid<Vec3> &gradient, RandomStream &rand)
   {
-    static RandomStream rand(3123984);
     const double jlen = 0.1;
 
     if (part.isActive(idx)) {
@@ -2321,6 +2330,11 @@ template<class S> struct KnProjectParticles : public KernelBase {
     return gradient;
   }
   typedef Grid<Vec3> type1;
+  inline RandomStream &getArg2()
+  {
+    return rand;
+  }
+  typedef RandomStream type2;
   void runMessage()
   {
     debMsg("Executing kernel KnProjectParticles ", 3);
@@ -2332,15 +2346,17 @@ template<class S> struct KnProjectParticles : public KernelBase {
   {
     const IndexInt _sz = size;
     for (IndexInt i = 0; i < _sz; i++)
-      op(i, part, gradient);
+      op(i, part, gradient, rand);
   }
   ParticleSystem<S> &part;
   Grid<Vec3> &gradient;
+  RandomStream &rand;
 };
 
 template<class S> void ParticleSystem<S>::projectOutside(Grid<Vec3> &gradient)
 {
-  KnProjectParticles<S>(*this, gradient);
+  RandomStream rand(globalSeed);
+  KnProjectParticles<S>(*this, gradient, rand);
 }
 
 template<class S> struct KnProjectOutOfBnd : public KernelBase {
@@ -2531,14 +2547,14 @@ template<class S> void ParticleSystem<S>::insertBufferedParticles()
 
   int insertFlag;
   Vec3 insertPos;
-  static RandomStream mRand(9832);
+  RandomStream rand(globalSeed);
   for (IndexInt i = 0; i < numNewParts; ++i) {
 
     // get random index in newBuffer vector
     // we are inserting particles randomly so that they are sampled uniformly in the fluid region
     // otherwise, regions of fluid can remain completely empty once mData.size() == maxParticles is
     // reached.
-    int randIndex = floor(mRand.getReal() * mNewBufferPos.size());
+    int randIndex = floor(rand.getReal() * mNewBufferPos.size());
 
     // get elements from new buffers with random index
     std::swap(mNewBufferPos[randIndex], mNewBufferPos.back());
diff --git a/extern/mantaflow/preprocessed/plugin/flip.cpp b/extern/mantaflow/preprocessed/plugin/flip.cpp
index 1acdac1c094..eaa1ebe45d5 100644
--- a/extern/mantaflow/preprocessed/plugin/flip.cpp
+++ b/extern/mantaflow/preprocessed/plugin/flip.cpp
@@ -41,7 +41,7 @@ void sampleFlagsWithParticles(const FlagGrid &flags,
   const bool is3D = flags.is3D();
   const Real jlen = randomness / discretization;
   const Vec3 disp(1.0 / discretization, 1.0 / discretization, 1.0 / discretization);
-  RandomStream mRand(9832);
+  RandomStream rand(parts.getSeed());
 
   FOR_IJK_BND(flags, 0)
   {
@@ -53,7 +53,7 @@ void sampleFlagsWithParticles(const FlagGrid &flags,
         for (int dj = 0; dj < discretization; dj++)
           for (int di = 0; di < discretization; di++) {
             Vec3 subpos = pos + disp * Vec3(0.5 + di, 0.5 + dj, 0.5 + dk);
-            subpos += jlen * (Vec3(1, 1, 1) - 2.0 * mRand.getVec3());
+            subpos += jlen * (Vec3(1, 1, 1) - 2.0 * rand.getVec3());
             if (!is3D)
               subpos[2] = 0.5;
             parts.addBuffered(subpos);
@@ -113,7 +113,7 @@ void sampleLevelsetWithParticles(const LevelsetGrid &phi,
   const bool is3D = phi.is3D();
   const Real jlen = randomness / discretization;
   const Vec3 disp(1.0 / discretization, 1.0 / discretization, 1.0 / discretization);
-  RandomStream mRand(9832);
+  RandomStream rand(parts.getSeed());
 
   if (reset) {
     parts.clear();
@@ -132,7 +132,7 @@ void sampleLevelsetWithParticles(const LevelsetGrid &phi,
         for (int dj = 0; dj < discretization; dj++)
           for (int di = 0; di < discretization; di++) {
             Vec3 subpos = pos + disp * Vec3(0.5 + di, 0.5 + dj, 0.5 + dk);
-            subpos += jlen * (Vec3(1, 1, 1) - 2.0 * mRand.getVec3());
+            subpos += jlen * (Vec3(1, 1, 1) - 2.0 * rand.getVec3());
             if (!is3D)
               subpos[2] = 0.5;
             if (phi.getInterpolated(subpos) > 0.)
@@ -205,7 +205,7 @@ void sampleShapeWithParticles(const Shape &shape,
   const bool is3D = flags.is3D();
   const Real jlen = randomness / discretization;
   const Vec3 disp(1.0 / discretization, 1.0 / discretization, 1.0 / discretization);
-  RandomStream mRand(9832);
+  RandomStream rand(parts.getSeed());
 
   if (reset) {
     parts.clear();
@@ -223,7 +223,7 @@ void sampleShapeWithParticles(const Shape &shape,
       for (int dj = 0; dj < discretization; dj++)
         for (int di = 0; di < discretization; di++) {
           Vec3 subpos = pos + disp * Vec3(0.5 + di, 0.5 + dj, 0.5 + dk);
-          subpos += jlen * (Vec3(1, 1, 1) - 2.0 * mRand.getVec3());
+          subpos += jlen * (Vec3(1, 1, 1) - 2.0 * rand.getVec3());
           if (!is3D)
             subpos[2] = 0.5;
           if (exclude && exclude->getInterpolated(subpos) <= 0.)
@@ -576,7 +576,7 @@ void adjustNumber(BasicParticleSystem &parts,
   }
 
   // seed new particles
-  RandomStream mRand(9832);
+  RandomStream rand(parts.getSeed());
   FOR_IJK(tmp)
   {
     int cnt = tmp(i, j, k);
@@ -593,7 +593,7 @@ void adjustNumber(BasicParticleSystem &parts,
 
     if (flags.isFluid(i, j, k) && cnt < minParticles) {
       for (int m = cnt; m < minParticles; m++) {
-        Vec3 pos = Vec3(i, j, k) + mRand.getVec3();
+        Vec3 pos = Vec3(i, j, k) + rand.getVec3();
         // Vec3 pos (i + 0.5, j + 0.5, k + 0.5); // cell center
         parts.addBuffered(pos);
       }
diff --git a/extern/mantaflow/preprocessed/plugin/secondaryparticles.cpp b/extern/mantaflow/preprocessed/plugin/secondaryparticles.cpp
index 82afe2d7ab2..8ebc239e8fc 100644
--- a/extern/mantaflow/preprocessed/plugin/secondaryparticles.cpp
+++ b/extern/mantaflow/preprocessed/plugin/secondaryparticles.cpp
@@ -479,7 +479,8 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase {
                                               const Real k_ta,
                                               const Real k_wc,
                                               const Real dt,
-                                              const int itype = FlagGrid::TypeFluid)
+                                              const int itype,
+                                              RandomStream &rand)
       : KernelBase(&flags, 0),
         flags(flags),
         v(v),
@@ -497,7 +498,8 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase {
         k_ta(k_ta),
         k_wc(k_wc),
         dt(dt),
-        itype(itype)
+        itype(itype),
+        rand(rand)
   {
     runMessage();
     run();
@@ -521,13 +523,13 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase {
                  const Real k_ta,
                  const Real k_wc,
                  const Real dt,
-                 const int itype = FlagGrid::TypeFluid)
+                 const int itype,
+                 RandomStream &rand)
   {
 
     if (!(flags(i, j, k) & itype))
       return;
 
-    static RandomStream mRand(9832);
     Real radius =
         0.25;  // diameter=0.5 => sampling with two cylinders in each dimension since cell size=1
     for (Real x = i - radius; x <= i + radius; x += 2 * radius) {
@@ -549,9 +551,9 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase {
               cross(e1, dir));  // perpendicular to dir and e1, so e1 and e1 create reference plane
 
           for (int di = 0; di < n; di++) {
-            const Real r = radius * sqrt(mRand.getReal());        // distance to cylinder axis
-            const Real theta = mRand.getReal() * Real(2) * M_PI;  // azimuth
-            const Real h = mRand.getReal() * norm(dt * vi);       // distance to reference plane
+            const Real r = radius * sqrt(rand.getReal());        // distance to cylinder axis
+            const Real theta = rand.getReal() * Real(2) * M_PI;  // azimuth
+            const Real h = rand.getReal() * norm(dt * vi);       // distance to reference plane
             Vec3 xd = xi + r * cos(theta) * e1 + r * sin(theta) * e2 + h * getNormalized(vi);
             if (!flags.is3D())
               xd.z = 0;
@@ -561,7 +563,7 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase {
                                       vi;  // init velocity of new particle
             Real temp = (KE + TA + WC) / 3;
             l_sec[l_sec.size() - 1] = ((lMax - lMin) * temp) + lMin +
-                                      mRand.getReal() * 0.1;  // init lifetime of new particle
+                                      rand.getReal() * 0.1;  // init lifetime of new particle
 
             // init type of new particle
             if (neighborRatio(i, j, k) < c_s) {
@@ -663,6 +665,11 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase {
     return itype;
   }
   typedef int type16;
+  inline RandomStream &getArg17()
+  {
+    return rand;
+  }
+  typedef RandomStream type17;
   void runMessage()
   {
     debMsg("Executing kernel knFlipSampleSecondaryParticlesMoreCylinders ", 3);
@@ -696,7 +703,8 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase {
              k_ta,
              k_wc,
              dt,
-             itype);
+             itype,
+             rand);
   }
   const FlagGrid &flags;
   const MACGrid &v;
@@ -715,6 +723,7 @@ struct knFlipSampleSecondaryParticlesMoreCylinders : public KernelBase {
   const Real k_wc;
   const Real dt;
   const int itype;
+  RandomStream &rand;
 };
 
 // adds secondary particles to &pts_sec for every fluid cell in &flags according to the potential
@@ -738,7 +747,8 @@ struct knFlipSampleSecondaryParticles : public KernelBase {
                                  const Real k_ta,
                                  const Real k_wc,
                                  const Real dt,
-                                 const int itype = FlagGrid::TypeFluid)
+                                 const int itype,
+                                 RandomStream &rand)
       : KernelBase(&flags, 0),
         flags(flags),
         v(v),
@@ -756,7 +766,8 @@ struct knFlipSampleSecondaryParticles : public KernelBase {
         k_ta(k_ta),
         k_wc(k_wc),
         dt(dt),
-        itype(itype)
+        itype(itype),
+        rand(rand)
   {
     runMessage();
     run();
@@ -780,7 +791,8 @@ struct knFlipSampleSecondaryParticles : public KernelBase {
                  const Real k_ta,
                  const Real k_wc,
                  const Real dt,
-                 const int itype = FlagGrid::TypeFluid)
+                 const int itype,
+                 RandomStream &rand)
   {
 
     if (!(flags(i, j, k) & itype))
@@ -793,9 +805,8 @@ struct knFlipSampleSecondaryParticles : public KernelBase {
     const int n = KE * (k_ta * TA + k_wc * WC) * dt;  // number of secondary particles
     if (n == 0)
       return;
-    static RandomStream mRand(9832);
 
-    Vec3 xi = Vec3(i, j, k) + mRand.getVec3();  // randomized offset uniform in cell
+    Vec3 xi = Vec3(i, j, k) + rand.getVec3();  // randomized offset uniform in cell
     Vec3 vi = v.getInterpolated(xi);
     Vec3 dir = dt * vi;                               // direction of movement of current particle
     Vec3 e1 = getNormalized(Vec3(dir.z, 0, -dir.x));  // perpendicular to dir
@@ -803,9 +814,9 @@ struct knFlipSampleSecondaryParticles : public KernelBase {
         cross(e1, dir));  // perpendicular to dir and e1, so e1 and e1 create reference plane
 
     for (int di = 0; di < n; di++) {
-      const Real r = Real(0.5) * sqrt(mRand.getReal());     // distance to cylinder axis
-      const Real theta = mRand.getReal() * Real(2) * M_PI;  // azimuth
-      const Real h = mRand.getReal() * norm(dt * vi);       // distance to reference plane
+      const Real r = Real(0.5) * sqrt(rand.getReal());     // distance to cylinder axis
+      const Real theta = rand.getReal() * Real(2) * M_PI;  // azimuth
+      const Real h = rand.getReal() * norm(dt * vi);       // distance to reference plane
       Vec3 xd = xi + r * cos(theta) * e1 + r * sin(theta) * e2 + h * getNormalized(vi);
       if (!flags.is3D())
         xd.z = 0;
@@ -815,7 +826,7 @@ struct knFlipSampleSecondaryParticles : public KernelBase {
                                 vi;  // init velocity of new particle
       Real temp = (KE + TA + WC) / 3;
       l_sec[l_sec.size() - 1] = ((lMax - lMin) * temp) + lMin +
-                                mRand.getReal() * 0.1;  // init lifetime of new particle
+                                rand.getReal() * 0.1;  // init lifetime of new particle
 
       // init type of new particle
       if (neighborRatio(i, j, k) < c_s) {
@@ -914,6 +925,11 @@ struct knFlipSampleSecondaryParticles : public KernelBase {
     return itype;
   }
   typedef int type16;
+  inline RandomStream &getArg17()
+  {
+    return rand;
+  }
+  typedef RandomStream type17;
   void runMessage()
   {
     debMsg("Executing kernel knFlipSampleSecondaryParticles ", 3);
@@ -947,7 +963,8 @@ struct knFlipSampleSecondaryParticles : public KernelBase {
              k_ta,
              k_wc,
              dt,
-             itype);
+             itype,
+             rand);
   }
   const FlagGrid &flags;
   const MACGrid &v;
@@ -966,6 +983,7 @@ struct knFlipSampleSecondaryParticles : public KernelBase {
   const Real k_wc;
   const Real dt;
   const int itype;
+  RandomStream &rand;
 };
 
 void flipSampleSecondaryParticles(const std::string mode,
@@ -992,6 +1010,9 @@ void flipSampleSecondaryParticles(const std::string mode,
   if (dt <= 0)
     timestep = flags.getParent()->getDt();
 
+  /* Every particle needs to get a different random offset. */
+  RandomStream rand(pts_sec.getSeed());
+
   if (mode == "single") {
     knFlipSampleSecondaryParticles(flags,
                                    v,
@@ -1009,7 +1030,8 @@ void flipSampleSecondaryParticles(const std::string mode,
                                    k_ta,
                                    k_wc,
                                    timestep,
-                                   itype);
+                                   itype,
+                                   rand);
   }
   else if (mode == "multiple") {
     knFlipSampleSecondaryParticlesMoreCylinders(flags,
@@ -1028,7 +1050,8 @@ void flipSampleSecondaryParticles(const std::string mode,
                                                 k_ta,
                                                 k_wc,
                                                 timestep,
-                                                itype);
+                                                itype,
+                                                rand);
   }
   else {
     throw std::invalid_argument("Unknown mode: use \"single\" or \"multiple\" instead!");