Welcome to mirror list, hosted at ThFree Co, Russian Federation.

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
path: root/extern/mantaflow/preprocessed/plugin/waves.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'extern/mantaflow/preprocessed/plugin/waves.cpp')
-rw-r--r--extern/mantaflow/preprocessed/plugin/waves.cpp483
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
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-06 16:11:58 +0300