Mantaflow [Part 1]: Added preprocessed Mantaflow source files

Includes preprocessed Mantaflow source files for both OpenMP and TBB (if OpenMP is not present, TBB files will be used instead).

These files come directly from the Mantaflow repository. Future updates to the core fluid solver will take place by updating the files.

Reviewed By: sergey, mont29

Maniphest Tasks: T59995

Differential Revision: https://developer.blender.org/D3850

1 files changed, 483 insertions, 0 deletions

diff --git a/extern/mantaflow/preprocessed/plugin/waves.cpp b/extern/mantaflow/preprocessed/plugin/waves.cpp
new file mode 100644
index 00000000000..7745dce4711
--- /dev/null
+++ b/extern/mantaflow/preprocessed/plugin/waves.cpp
@@ -0,0 +1,483 @@
+
+
+// DO NOT EDIT !
+// This file is generated using the MantaFlow preprocessor (prep generate).
+
+/******************************************************************************
+ *
+ * MantaFlow fluid solver framework
+ * Copyright 2011 Tobias Pfaff, 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
+ *
+ * Wave equation
+ *
+ ******************************************************************************/
+
+#include "levelset.h"
+#include "commonkernels.h"
+#include "particle.h"
+#include "conjugategrad.h"
+#include <cmath>
+
+using namespace std;
+
+namespace Manta {
+
+/******************************************************************************
+ *
+ * explicit integration
+ *
+ ******************************************************************************/
+
+struct knCalcSecDeriv2d : public KernelBase {
+  knCalcSecDeriv2d(const Grid<Real> &v, Grid<Real> &ret) : KernelBase(&v, 1), v(v), ret(ret)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(int i, int j, int k, const Grid<Real> &v, Grid<Real> &ret) const
+  {
+    ret(i, j, k) = (-4. * v(i, j, k) + v(i - 1, j, k) + v(i + 1, j, k) + v(i, j - 1, k) +
+                    v(i, j + 1, k));
+  }
+  inline const Grid<Real> &getArg0()
+  {
+    return v;
+  }
+  typedef Grid<Real> type0;
+  inline Grid<Real> &getArg1()
+  {
+    return ret;
+  }
+  typedef Grid<Real> type1;
+  void runMessage()
+  {
+    debMsg("Executing kernel knCalcSecDeriv2d ", 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 = 1; j < _maxY; j++)
+          for (int i = 1; i < _maxX; i++)
+            op(i, j, k, v, ret);
+    }
+    else {
+      const int k = 0;
+      for (int j = __r.begin(); j != (int)__r.end(); j++)
+        for (int i = 1; i < _maxX; i++)
+          op(i, j, k, v, ret);
+    }
+  }
+  void run()
+  {
+    if (maxZ > 1)
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
+    else
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
+  }
+  const Grid<Real> &v;
+  Grid<Real> &ret;
+};
+;
+
+//! calculate a second derivative for the wave equation
+void calcSecDeriv2d(const Grid<Real> &v, Grid<Real> &curv)
+{
+  knCalcSecDeriv2d(v, curv);
+}
+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, "calcSecDeriv2d", !noTiming);
+    PyObject *_retval = 0;
+    {
+      ArgLocker _lock;
+      const Grid<Real> &v = *_args.getPtr<Grid<Real>>("v", 0, &_lock);
+      Grid<Real> &curv = *_args.getPtr<Grid<Real>>("curv", 1, &_lock);
+      _retval = getPyNone();
+      calcSecDeriv2d(v, curv);
+      _args.check();
+    }
+    pbFinalizePlugin(parent, "calcSecDeriv2d", !noTiming);
+    return _retval;
+  }
+  catch (std::exception &e) {
+    pbSetError("calcSecDeriv2d", e.what());
+    return 0;
+  }
+}
+static const Pb::Register _RP_calcSecDeriv2d("", "calcSecDeriv2d", _W_0);
+extern "C" {
+void PbRegister_calcSecDeriv2d()
+{
+  KEEP_UNUSED(_RP_calcSecDeriv2d);
+}
+}
+
+// mass conservation
+
+struct knTotalSum : public KernelBase {
+  knTotalSum(Grid<Real> &h) : KernelBase(&h, 1), h(h), sum(0)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(int i, int j, int k, Grid<Real> &h, double &sum)
+  {
+    sum += h(i, j, k);
+  }
+  inline operator double()
+  {
+    return sum;
+  }
+  inline double &getRet()
+  {
+    return sum;
+  }
+  inline Grid<Real> &getArg0()
+  {
+    return h;
+  }
+  typedef Grid<Real> type0;
+  void runMessage()
+  {
+    debMsg("Executing kernel knTotalSum ", 3);
+    debMsg("Kernel range"
+               << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ",
+           4);
+  };
+  void operator()(const tbb::blocked_range<IndexInt> &__r)
+  {
+    const int _maxX = maxX;
+    const int _maxY = maxY;
+    if (maxZ > 1) {
+      for (int k = __r.begin(); k != (int)__r.end(); k++)
+        for (int j = 1; j < _maxY; j++)
+          for (int i = 1; i < _maxX; i++)
+            op(i, j, k, h, sum);
+    }
+    else {
+      const int k = 0;
+      for (int j = __r.begin(); j != (int)__r.end(); j++)
+        for (int i = 1; i < _maxX; i++)
+          op(i, j, k, h, sum);
+    }
+  }
+  void run()
+  {
+    if (maxZ > 1)
+      tbb::parallel_reduce(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
+    else
+      tbb::parallel_reduce(tbb::blocked_range<IndexInt>(1, maxY), *this);
+  }
+  knTotalSum(knTotalSum &o, tbb::split) : KernelBase(o), h(o.h), sum(0)
+  {
+  }
+  void join(const knTotalSum &o)
+  {
+    sum += o.sum;
+  }
+  Grid<Real> &h;
+  double sum;
+};
+
+//! calculate the sum of all values in a grid (for wave equation solves)
+Real totalSum(Grid<Real> &height)
+{
+  knTotalSum ts(height);
+  return ts.sum;
+}
+static PyObject *_W_1(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, "totalSum", !noTiming);
+    PyObject *_retval = 0;
+    {
+      ArgLocker _lock;
+      Grid<Real> &height = *_args.getPtr<Grid<Real>>("height", 0, &_lock);
+      _retval = toPy(totalSum(height));
+      _args.check();
+    }
+    pbFinalizePlugin(parent, "totalSum", !noTiming);
+    return _retval;
+  }
+  catch (std::exception &e) {
+    pbSetError("totalSum", e.what());
+    return 0;
+  }
+}
+static const Pb::Register _RP_totalSum("", "totalSum", _W_1);
+extern "C" {
+void PbRegister_totalSum()
+{
+  KEEP_UNUSED(_RP_totalSum);
+}
+}
+
+//! normalize all values in a grid (for wave equation solves)
+void normalizeSumTo(Grid<Real> &height, Real target)
+{
+  knTotalSum ts(height);
+  Real factor = target / ts.sum;
+  height.multConst(factor);
+}
+static PyObject *_W_2(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, "normalizeSumTo", !noTiming);
+    PyObject *_retval = 0;
+    {
+      ArgLocker _lock;
+      Grid<Real> &height = *_args.getPtr<Grid<Real>>("height", 0, &_lock);
+      Real target = _args.get<Real>("target", 1, &_lock);
+      _retval = getPyNone();
+      normalizeSumTo(height, target);
+      _args.check();
+    }
+    pbFinalizePlugin(parent, "normalizeSumTo", !noTiming);
+    return _retval;
+  }
+  catch (std::exception &e) {
+    pbSetError("normalizeSumTo", e.what());
+    return 0;
+  }
+}
+static const Pb::Register _RP_normalizeSumTo("", "normalizeSumTo", _W_2);
+extern "C" {
+void PbRegister_normalizeSumTo()
+{
+  KEEP_UNUSED(_RP_normalizeSumTo);
+}
+}
+
+/******************************************************************************
+ *
+ * implicit time integration
+ *
+ ******************************************************************************/
+
+//! Kernel: Construct the right-hand side of the poisson equation
+
+struct MakeRhsWE : public KernelBase {
+  MakeRhsWE(const FlagGrid &flags,
+            Grid<Real> &rhs,
+            const Grid<Real> &ut,
+            const Grid<Real> &utm1,
+            Real s,
+            bool crankNic = false)
+      : KernelBase(&flags, 1), flags(flags), rhs(rhs), ut(ut), utm1(utm1), s(s), crankNic(crankNic)
+  {
+    runMessage();
+    run();
+  }
+  inline void op(int i,
+                 int j,
+                 int k,
+                 const FlagGrid &flags,
+                 Grid<Real> &rhs,
+                 const Grid<Real> &ut,
+                 const Grid<Real> &utm1,
+                 Real s,
+                 bool crankNic = false) const
+  {
+    rhs(i, j, k) = (2. * ut(i, j, k) - utm1(i, j, k));
+    if (crankNic) {
+      rhs(i, j, k) += s * (-4. * ut(i, j, k) + 1. * ut(i - 1, j, k) + 1. * ut(i + 1, j, k) +
+                           1. * ut(i, j - 1, k) + 1. * ut(i, j + 1, k));
+    }
+  }
+  inline const FlagGrid &getArg0()
+  {
+    return flags;
+  }
+  typedef FlagGrid type0;
+  inline Grid<Real> &getArg1()
+  {
+    return rhs;
+  }
+  typedef Grid<Real> type1;
+  inline const Grid<Real> &getArg2()
+  {
+    return ut;
+  }
+  typedef Grid<Real> type2;
+  inline const Grid<Real> &getArg3()
+  {
+    return utm1;
+  }
+  typedef Grid<Real> type3;
+  inline Real &getArg4()
+  {
+    return s;
+  }
+  typedef Real type4;
+  inline bool &getArg5()
+  {
+    return crankNic;
+  }
+  typedef bool type5;
+  void runMessage()
+  {
+    debMsg("Executing kernel MakeRhsWE ", 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 = 1; j < _maxY; j++)
+          for (int i = 1; i < _maxX; i++)
+            op(i, j, k, flags, rhs, ut, utm1, s, crankNic);
+    }
+    else {
+      const int k = 0;
+      for (int j = __r.begin(); j != (int)__r.end(); j++)
+        for (int i = 1; i < _maxX; i++)
+          op(i, j, k, flags, rhs, ut, utm1, s, crankNic);
+    }
+  }
+  void run()
+  {
+    if (maxZ > 1)
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(minZ, maxZ), *this);
+    else
+      tbb::parallel_for(tbb::blocked_range<IndexInt>(1, maxY), *this);
+  }
+  const FlagGrid &flags;
+  Grid<Real> &rhs;
+  const Grid<Real> &ut;
+  const Grid<Real> &utm1;
+  Real s;
+  bool crankNic;
+};
+
+//! do a CG solve for the wave equation (note, out grid only there for debugging... could be
+//! removed)
+
+void cgSolveWE(const FlagGrid &flags,
+               Grid<Real> &ut,
+               Grid<Real> &utm1,
+               Grid<Real> &out,
+               bool crankNic = false,
+               Real cSqr = 0.25,
+               Real cgMaxIterFac = 1.5,
+               Real cgAccuracy = 1e-5)
+{
+  // reserve temp grids
+  FluidSolver *parent = flags.getParent();
+  Grid<Real> rhs(parent);
+  Grid<Real> residual(parent);
+  Grid<Real> search(parent);
+  Grid<Real> A0(parent);
+  Grid<Real> Ai(parent);
+  Grid<Real> Aj(parent);
+  Grid<Real> Ak(parent);
+  Grid<Real> tmp(parent);
+  // solution...
+  out.clear();
+
+  // setup matrix and boundaries
+  MakeLaplaceMatrix(flags, A0, Ai, Aj, Ak);
+  Real dt = parent->getDt();
+  Real s = dt * dt * cSqr * 0.5;
+  FOR_IJK(flags)
+  {
+    Ai(i, j, k) *= s;
+    Aj(i, j, k) *= s;
+    Ak(i, j, k) *= s;
+    A0(i, j, k) *= s;
+    A0(i, j, k) += 1.;
+  }
+
+  // compute divergence and init right hand side
+  rhs.clear();
+  // h=dt
+  // rhs:   = 2 ut - ut-1
+  // A:    (h2 c2/ dx)=s   ,  (1+4s)uij + s ui-1j + ...
+  // Cr.Nic.
+  // rhs:  cr nic = 2 ut - ut-1 + h^2c^2/2 b
+  // A:    (h2 c2/2 dx)=s   ,  (1+4s)uij + s ui-1j + ...
+  MakeRhsWE kernMakeRhs(flags, rhs, ut, utm1, s, crankNic);
+
+  const int maxIter = (int)(cgMaxIterFac * flags.getSize().max()) * (flags.is3D() ? 1 : 4);
+  GridCgInterface *gcg;
+  if (flags.is3D())
+    gcg = new GridCg<ApplyMatrix>(out, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+  else
+    gcg = new GridCg<ApplyMatrix2D>(out, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak);
+
+  gcg->setAccuracy(cgAccuracy);
+
+  // no preconditioning for now...
+  for (int iter = 0; iter < maxIter; iter++) {
+    if (!gcg->iterate())
+      iter = maxIter;
+  }
+  debMsg("cgSolveWaveEq iterations:" << gcg->getIterations() << ", res:" << gcg->getSigma(), 1);
+
+  utm1.swap(ut);
+  ut.copyFrom(out);
+
+  delete gcg;
+}
+static PyObject *_W_3(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, "cgSolveWE", !noTiming);
+    PyObject *_retval = 0;
+    {
+      ArgLocker _lock;
+      const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 0, &_lock);
+      Grid<Real> &ut = *_args.getPtr<Grid<Real>>("ut", 1, &_lock);
+      Grid<Real> &utm1 = *_args.getPtr<Grid<Real>>("utm1", 2, &_lock);
+      Grid<Real> &out = *_args.getPtr<Grid<Real>>("out", 3, &_lock);
+      bool crankNic = _args.getOpt<bool>("crankNic", 4, false, &_lock);
+      Real cSqr = _args.getOpt<Real>("cSqr", 5, 0.25, &_lock);
+      Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 6, 1.5, &_lock);
+      Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 7, 1e-5, &_lock);
+      _retval = getPyNone();
+      cgSolveWE(flags, ut, utm1, out, crankNic, cSqr, cgMaxIterFac, cgAccuracy);
+      _args.check();
+    }
+    pbFinalizePlugin(parent, "cgSolveWE", !noTiming);
+    return _retval;
+  }
+  catch (std::exception &e) {
+    pbSetError("cgSolveWE", e.what());
+    return 0;
+  }
+}
+static const Pb::Register _RP_cgSolveWE("", "cgSolveWE", _W_3);
+extern "C" {
+void PbRegister_cgSolveWE()
+{
+  KEEP_UNUSED(_RP_cgSolveWE);
+}
+}
+
+}  // namespace Manta