diff options
Diffstat (limited to 'extern/mantaflow/preprocessed/plugin/pressure.cpp')
| -rw-r--r-- | extern/mantaflow/preprocessed/plugin/pressure.cpp | 1511 |
1 files changed, 1511 insertions, 0 deletions
diff --git a/extern/mantaflow/preprocessed/plugin/pressure.cpp b/extern/mantaflow/preprocessed/plugin/pressure.cpp new file mode 100644 index 00000000000..7def2669e36 --- /dev/null +++ b/extern/mantaflow/preprocessed/plugin/pressure.cpp @@ -0,0 +1,1511 @@ + + +// 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 + * + * Plugins for pressure correction: solve_pressure, and ghost fluid helpers + * + ******************************************************************************/ +#include "vectorbase.h" +#include "kernel.h" +#include "conjugategrad.h" +#include "multigrid.h" + +using namespace std; +namespace Manta { + +//! Preconditioner for CG solver +// - None: Use standard CG +// - MIC: Modified incomplete Cholesky preconditioner +// - MGDynamic: Multigrid preconditioner, rebuilt for each solve +// - MGStatic: Multigrid preconditioner, built only once (faster than +// MGDynamic, but works only if Poisson equation does not change) +enum Preconditioner { PcNone = 0, PcMIC = 1, PcMGDynamic = 2, PcMGStatic = 3 }; + +inline static Real surfTensHelper(const IndexInt idx, + const int offset, + const Grid<Real> &phi, + const Grid<Real> &curv, + const Real surfTens, + const Real gfClamp); + +//! Kernel: Construct the right-hand side of the poisson equation + +struct MakeRhs : public KernelBase { + MakeRhs(const FlagGrid &flags, + Grid<Real> &rhs, + const MACGrid &vel, + const Grid<Real> *perCellCorr, + const MACGrid *fractions, + const MACGrid *obvel, + const Grid<Real> *phi, + const Grid<Real> *curv, + const Real surfTens, + const Real gfClamp) + : KernelBase(&flags, 1), + flags(flags), + rhs(rhs), + vel(vel), + perCellCorr(perCellCorr), + fractions(fractions), + obvel(obvel), + phi(phi), + curv(curv), + surfTens(surfTens), + gfClamp(gfClamp), + cnt(0), + sum(0) + { + runMessage(); + run(); + } + inline void op(int i, + int j, + int k, + const FlagGrid &flags, + Grid<Real> &rhs, + const MACGrid &vel, + const Grid<Real> *perCellCorr, + const MACGrid *fractions, + const MACGrid *obvel, + const Grid<Real> *phi, + const Grid<Real> *curv, + const Real surfTens, + const Real gfClamp, + int &cnt, + double &sum) + { + if (!flags.isFluid(i, j, k)) { + rhs(i, j, k) = 0; + return; + } + + // compute negative divergence + // no flag checks: assumes vel at obstacle interfaces is set to zero + Real set(0); + if (!fractions) { + set = vel(i, j, k).x - vel(i + 1, j, k).x + vel(i, j, k).y - vel(i, j + 1, k).y; + if (vel.is3D()) + set += vel(i, j, k).z - vel(i, j, k + 1).z; + } + else { + set = (*fractions)(i, j, k).x * vel(i, j, k).x - + (*fractions)(i + 1, j, k).x * vel(i + 1, j, k).x + + (*fractions)(i, j, k).y * vel(i, j, k).y - + (*fractions)(i, j + 1, k).y * vel(i, j + 1, k).y; + if (vel.is3D()) + set += (*fractions)(i, j, k).z * vel(i, j, k).z - + (*fractions)(i, j, k + 1).z * vel(i, j, k + 1).z; + + // compute divergence from obstacle by using obstacle velocity (optional) + if (obvel) { + set += (1 - (*fractions)(i, j, k).x) * (*obvel)(i, j, k).x - + (1 - (*fractions)(i + 1, j, k).x) * (*obvel)(i + 1, j, k).x + + (1 - (*fractions)(i, j, k).y) * (*obvel)(i, j, k).y - + (1 - (*fractions)(i, j + 1, k).y) * (*obvel)(i, j + 1, k).y; + if (obvel->is3D()) + set += (1 - (*fractions)(i, j, k).z) * (*obvel)(i, j, k).z - + (1 - (*fractions)(i, j, k + 1).z) * (*obvel)(i, j, k + 1).z; + } + } + + // compute surface tension effect (optional) + if (phi && curv) { + const IndexInt idx = flags.index(i, j, k); + const int X = flags.getStrideX(), Y = flags.getStrideY(), Z = flags.getStrideZ(); + if (flags.isEmpty(i - 1, j, k)) + set += surfTensHelper(idx, -X, *phi, *curv, surfTens, gfClamp); + if (flags.isEmpty(i + 1, j, k)) + set += surfTensHelper(idx, +X, *phi, *curv, surfTens, gfClamp); + if (flags.isEmpty(i, j - 1, k)) + set += surfTensHelper(idx, -Y, *phi, *curv, surfTens, gfClamp); + if (flags.isEmpty(i, j + 1, k)) + set += surfTensHelper(idx, +Y, *phi, *curv, surfTens, gfClamp); + if (vel.is3D()) { + if (flags.isEmpty(i, j, k - 1)) + set += surfTensHelper(idx, -Z, *phi, *curv, surfTens, gfClamp); + if (flags.isEmpty(i, j, k + 1)) + set += surfTensHelper(idx, +Z, *phi, *curv, surfTens, gfClamp); + } + } + + // per cell divergence correction (optional) + if (perCellCorr) + set += perCellCorr->get(i, j, k); + + // obtain sum, cell count + sum += set; + cnt++; + + rhs(i, j, k) = set; + } + inline const FlagGrid &getArg0() + { + return flags; + } + typedef FlagGrid type0; + inline Grid<Real> &getArg1() + { + return rhs; + } + typedef Grid<Real> type1; + inline const MACGrid &getArg2() + { + return vel; + } + typedef MACGrid type2; + inline const Grid<Real> *getArg3() + { + return perCellCorr; + } + typedef Grid<Real> type3; + inline const MACGrid *getArg4() + { + return fractions; + } + typedef MACGrid type4; + inline const MACGrid *getArg5() + { + return obvel; + } + typedef MACGrid type5; + inline const Grid<Real> *getArg6() + { + return phi; + } + typedef Grid<Real> type6; + inline const Grid<Real> *getArg7() + { + return curv; + } + typedef Grid<Real> type7; + inline const Real &getArg8() + { + return surfTens; + } + typedef Real type8; + inline const Real &getArg9() + { + return gfClamp; + } + typedef Real type9; + void runMessage() + { + debMsg("Executing kernel MakeRhs ", 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, + flags, + rhs, + vel, + perCellCorr, + fractions, + obvel, + phi, + curv, + surfTens, + gfClamp, + cnt, + 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, + flags, + rhs, + vel, + perCellCorr, + fractions, + obvel, + phi, + curv, + surfTens, + gfClamp, + cnt, + 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); + } + MakeRhs(MakeRhs &o, tbb::split) + : KernelBase(o), + flags(o.flags), + rhs(o.rhs), + vel(o.vel), + perCellCorr(o.perCellCorr), + fractions(o.fractions), + obvel(o.obvel), + phi(o.phi), + curv(o.curv), + surfTens(o.surfTens), + gfClamp(o.gfClamp), + cnt(0), + sum(0) + { + } + void join(const MakeRhs &o) + { + cnt += o.cnt; + sum += o.sum; + } + const FlagGrid &flags; + Grid<Real> &rhs; + const MACGrid &vel; + const Grid<Real> *perCellCorr; + const MACGrid *fractions; + const MACGrid *obvel; + const Grid<Real> *phi; + const Grid<Real> *curv; + const Real surfTens; + const Real gfClamp; + int cnt; + double sum; +}; + +//! Kernel: make velocity divergence free by subtracting pressure gradient + +struct knCorrectVelocity : public KernelBase { + knCorrectVelocity(const FlagGrid &flags, MACGrid &vel, const Grid<Real> &pressure) + : KernelBase(&flags, 1), flags(flags), vel(vel), pressure(pressure) + { + runMessage(); + run(); + } + inline void op( + int i, int j, int k, const FlagGrid &flags, MACGrid &vel, const Grid<Real> &pressure) const + { + const IndexInt idx = flags.index(i, j, k); + if (flags.isFluid(idx)) { + if (flags.isFluid(i - 1, j, k)) + vel[idx].x -= (pressure[idx] - pressure(i - 1, j, k)); + if (flags.isFluid(i, j - 1, k)) + vel[idx].y -= (pressure[idx] - pressure(i, j - 1, k)); + if (flags.is3D() && flags.isFluid(i, j, k - 1)) + vel[idx].z -= (pressure[idx] - pressure(i, j, k - 1)); + + if (flags.isEmpty(i - 1, j, k)) + vel[idx].x -= pressure[idx]; + if (flags.isEmpty(i, j - 1, k)) + vel[idx].y -= pressure[idx]; + if (flags.is3D() && flags.isEmpty(i, j, k - 1)) + vel[idx].z -= pressure[idx]; + } + else if (flags.isEmpty(idx) && + !flags.isOutflow(idx)) { // don't change velocities in outflow cells + if (flags.isFluid(i - 1, j, k)) + vel[idx].x += pressure(i - 1, j, k); + else + vel[idx].x = 0.f; + if (flags.isFluid(i, j - 1, k)) + vel[idx].y += pressure(i, j - 1, k); + else + vel[idx].y = 0.f; + if (flags.is3D()) { + if (flags.isFluid(i, j, k - 1)) + vel[idx].z += pressure(i, j, k - 1); + else + vel[idx].z = 0.f; + } + } + } + inline const FlagGrid &getArg0() + { + return flags; + } + typedef FlagGrid type0; + inline MACGrid &getArg1() + { + return vel; + } + typedef MACGrid type1; + inline const Grid<Real> &getArg2() + { + return pressure; + } + typedef Grid<Real> type2; + void runMessage() + { + debMsg("Executing kernel knCorrectVelocity ", 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, vel, pressure); + } + 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, vel, pressure); + } + } + 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; + MACGrid &vel; + const Grid<Real> &pressure; +}; + +// ***************************************************************************** +// Ghost fluid helpers + +// calculate fraction filled with liquid (note, assumes inside value is < outside!) +inline static Real thetaHelper(const Real inside, const Real outside) +{ + const Real denom = inside - outside; + if (denom > -1e-04) + return 0.5; // should always be neg, and large enough... + return std::max(Real(0), std::min(Real(1), inside / denom)); +} + +// calculate ghost fluid factor, cell at idx should be a fluid cell +inline static Real ghostFluidHelper(const IndexInt idx, + const int offset, + const Grid<Real> &phi, + const Real gfClamp) +{ + Real alpha = thetaHelper(phi[idx], phi[idx + offset]); + if (alpha < gfClamp) + return alpha = gfClamp; + return (1. - (1. / alpha)); +} + +inline static Real surfTensHelper(const IndexInt idx, + const int offset, + const Grid<Real> &phi, + const Grid<Real> &curv, + const Real surfTens, + const Real gfClamp) +{ + return surfTens * (curv[idx + offset] - ghostFluidHelper(idx, offset, phi, gfClamp) * curv[idx]); +} + +//! Kernel: Adapt A0 for ghost fluid + +struct ApplyGhostFluidDiagonal : public KernelBase { + ApplyGhostFluidDiagonal(Grid<Real> &A0, + const FlagGrid &flags, + const Grid<Real> &phi, + const Real gfClamp) + : KernelBase(&A0, 1), A0(A0), flags(flags), phi(phi), gfClamp(gfClamp) + { + runMessage(); + run(); + } + inline void op(int i, + int j, + int k, + Grid<Real> &A0, + const FlagGrid &flags, + const Grid<Real> &phi, + const Real gfClamp) const + { + const int X = flags.getStrideX(), Y = flags.getStrideY(), Z = flags.getStrideZ(); + const IndexInt idx = flags.index(i, j, k); + if (!flags.isFluid(idx)) + return; + + if (flags.isEmpty(i - 1, j, k)) + A0[idx] -= ghostFluidHelper(idx, -X, phi, gfClamp); + if (flags.isEmpty(i + 1, j, k)) + A0[idx] -= ghostFluidHelper(idx, +X, phi, gfClamp); + if (flags.isEmpty(i, j - 1, k)) + A0[idx] -= ghostFluidHelper(idx, -Y, phi, gfClamp); + if (flags.isEmpty(i, j + 1, k)) + A0[idx] -= ghostFluidHelper(idx, +Y, phi, gfClamp); + if (flags.is3D()) { + if (flags.isEmpty(i, j, k - 1)) + A0[idx] -= ghostFluidHelper(idx, -Z, phi, gfClamp); + if (flags.isEmpty(i, j, k + 1)) + A0[idx] -= ghostFluidHelper(idx, +Z, phi, gfClamp); + } + } + inline Grid<Real> &getArg0() + { + return A0; + } + typedef Grid<Real> type0; + inline const FlagGrid &getArg1() + { + return flags; + } + typedef FlagGrid type1; + inline const Grid<Real> &getArg2() + { + return phi; + } + typedef Grid<Real> type2; + inline const Real &getArg3() + { + return gfClamp; + } + typedef Real type3; + void runMessage() + { + debMsg("Executing kernel ApplyGhostFluidDiagonal ", 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, A0, flags, phi, gfClamp); + } + 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, A0, flags, phi, gfClamp); + } + } + 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); + } + Grid<Real> &A0; + const FlagGrid &flags; + const Grid<Real> φ + const Real gfClamp; +}; + +//! Kernel: Apply velocity update: ghost fluid contribution + +struct knCorrectVelocityGhostFluid : public KernelBase { + knCorrectVelocityGhostFluid(MACGrid &vel, + const FlagGrid &flags, + const Grid<Real> &pressure, + const Grid<Real> &phi, + Real gfClamp, + const Grid<Real> *curv, + const Real surfTens) + : KernelBase(&vel, 1), + vel(vel), + flags(flags), + pressure(pressure), + phi(phi), + gfClamp(gfClamp), + curv(curv), + surfTens(surfTens) + { + runMessage(); + run(); + } + inline void op(int i, + int j, + int k, + MACGrid &vel, + const FlagGrid &flags, + const Grid<Real> &pressure, + const Grid<Real> &phi, + Real gfClamp, + const Grid<Real> *curv, + const Real surfTens) const + { + const IndexInt X = flags.getStrideX(), Y = flags.getStrideY(), Z = flags.getStrideZ(); + const IndexInt idx = flags.index(i, j, k); + if (flags.isFluid(idx)) { + if (flags.isEmpty(i - 1, j, k)) + vel[idx][0] += pressure[idx] * ghostFluidHelper(idx, -X, phi, gfClamp); + if (flags.isEmpty(i, j - 1, k)) + vel[idx][1] += pressure[idx] * ghostFluidHelper(idx, -Y, phi, gfClamp); + if (flags.is3D() && flags.isEmpty(i, j, k - 1)) + vel[idx][2] += pressure[idx] * ghostFluidHelper(idx, -Z, phi, gfClamp); + } + else if (flags.isEmpty(idx) && + !flags.isOutflow(idx)) { // do not change velocities in outflow cells + if (flags.isFluid(i - 1, j, k)) + vel[idx][0] -= pressure(i - 1, j, k) * ghostFluidHelper(idx - X, +X, phi, gfClamp); + else + vel[idx].x = 0.f; + if (flags.isFluid(i, j - 1, k)) + vel[idx][1] -= pressure(i, j - 1, k) * ghostFluidHelper(idx - Y, +Y, phi, gfClamp); + else + vel[idx].y = 0.f; + if (flags.is3D()) { + if (flags.isFluid(i, j, k - 1)) + vel[idx][2] -= pressure(i, j, k - 1) * ghostFluidHelper(idx - Z, +Z, phi, gfClamp); + else + vel[idx].z = 0.f; + } + } + + if (curv) { + if (flags.isFluid(idx)) { + if (flags.isEmpty(i - 1, j, k)) + vel[idx].x += surfTensHelper(idx, -X, phi, *curv, surfTens, gfClamp); + if (flags.isEmpty(i, j - 1, k)) + vel[idx].y += surfTensHelper(idx, -Y, phi, *curv, surfTens, gfClamp); + if (flags.is3D() && flags.isEmpty(i, j, k - 1)) + vel[idx].z += surfTensHelper(idx, -Z, phi, *curv, surfTens, gfClamp); + } + else if (flags.isEmpty(idx) && + !flags.isOutflow(idx)) { // do not change velocities in outflow cells + vel[idx].x -= (flags.isFluid(i - 1, j, k)) ? + surfTensHelper(idx - X, +X, phi, *curv, surfTens, gfClamp) : + 0.f; + vel[idx].y -= (flags.isFluid(i, j - 1, k)) ? + surfTensHelper(idx - Y, +Y, phi, *curv, surfTens, gfClamp) : + 0.f; + if (flags.is3D()) + vel[idx].z -= (flags.isFluid(i, j, k - 1)) ? + surfTensHelper(idx - Z, +Z, phi, *curv, surfTens, gfClamp) : + 0.f; + } + } + } + inline MACGrid &getArg0() + { + return vel; + } + typedef MACGrid type0; + inline const FlagGrid &getArg1() + { + return flags; + } + typedef FlagGrid type1; + inline const Grid<Real> &getArg2() + { + return pressure; + } + typedef Grid<Real> type2; + inline const Grid<Real> &getArg3() + { + return phi; + } + typedef Grid<Real> type3; + inline Real &getArg4() + { + return gfClamp; + } + typedef Real type4; + inline const Grid<Real> *getArg5() + { + return curv; + } + typedef Grid<Real> type5; + inline const Real &getArg6() + { + return surfTens; + } + typedef Real type6; + void runMessage() + { + debMsg("Executing kernel knCorrectVelocityGhostFluid ", 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, vel, flags, pressure, phi, gfClamp, curv, surfTens); + } + 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, vel, flags, pressure, phi, gfClamp, curv, surfTens); + } + } + 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); + } + MACGrid &vel; + const FlagGrid &flags; + const Grid<Real> &pressure; + const Grid<Real> φ + Real gfClamp; + const Grid<Real> *curv; + const Real surfTens; +}; + +// improve behavior of clamping for large time steps: +inline static Real ghostFluidWasClamped(const IndexInt idx, + const int offset, + const Grid<Real> &phi, + const Real gfClamp) +{ + const Real alpha = thetaHelper(phi[idx], phi[idx + offset]); + if (alpha < gfClamp) + return true; + return false; +} + +struct knReplaceClampedGhostFluidVels : public KernelBase { + knReplaceClampedGhostFluidVels(MACGrid &vel, + const FlagGrid &flags, + const Grid<Real> &pressure, + const Grid<Real> &phi, + Real gfClamp) + : KernelBase(&vel, 1), vel(vel), flags(flags), pressure(pressure), phi(phi), gfClamp(gfClamp) + { + runMessage(); + run(); + } + inline void op(int i, + int j, + int k, + MACGrid &vel, + const FlagGrid &flags, + const Grid<Real> &pressure, + const Grid<Real> &phi, + Real gfClamp) const + { + const IndexInt idx = flags.index(i, j, k); + const IndexInt X = flags.getStrideX(), Y = flags.getStrideY(), Z = flags.getStrideZ(); + if (!flags.isEmpty(idx)) + return; + + if (flags.isFluid(i - 1, j, k) && ghostFluidWasClamped(idx - X, +X, phi, gfClamp)) + vel[idx][0] = vel[idx - X][0]; + if (flags.isFluid(i, j - 1, k) && ghostFluidWasClamped(idx - Y, +Y, phi, gfClamp)) + vel[idx][1] = vel[idx - Y][1]; + if (flags.is3D() && flags.isFluid(i, j, k - 1) && + ghostFluidWasClamped(idx - Z, +Z, phi, gfClamp)) + vel[idx][2] = vel[idx - Z][2]; + + if (flags.isFluid(i + 1, j, k) && ghostFluidWasClamped(idx + X, -X, phi, gfClamp)) + vel[idx][0] = vel[idx + X][0]; + if (flags.isFluid(i, j + 1, k) && ghostFluidWasClamped(idx + Y, -Y, phi, gfClamp)) + vel[idx][1] = vel[idx + Y][1]; + if (flags.is3D() && flags.isFluid(i, j, k + 1) && + ghostFluidWasClamped(idx + Z, -Z, phi, gfClamp)) + vel[idx][2] = vel[idx + Z][2]; + } + inline MACGrid &getArg0() + { + return vel; + } + typedef MACGrid type0; + inline const FlagGrid &getArg1() + { + return flags; + } + typedef FlagGrid type1; + inline const Grid<Real> &getArg2() + { + return pressure; + } + typedef Grid<Real> type2; + inline const Grid<Real> &getArg3() + { + return phi; + } + typedef Grid<Real> type3; + inline Real &getArg4() + { + return gfClamp; + } + typedef Real type4; + void runMessage() + { + debMsg("Executing kernel knReplaceClampedGhostFluidVels ", 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, vel, flags, pressure, phi, gfClamp); + } + 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, vel, flags, pressure, phi, gfClamp); + } + } + 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); + } + MACGrid &vel; + const FlagGrid &flags; + const Grid<Real> &pressure; + const Grid<Real> φ + Real gfClamp; +}; + +//! Kernel: Compute min value of Real grid + +struct CountEmptyCells : public KernelBase { + CountEmptyCells(const FlagGrid &flags) : KernelBase(&flags, 0), flags(flags), numEmpty(0) + { + runMessage(); + run(); + } + inline void op(IndexInt idx, const FlagGrid &flags, int &numEmpty) + { + if (flags.isEmpty(idx)) + numEmpty++; + } + inline operator int() + { + return numEmpty; + } + inline int &getRet() + { + return numEmpty; + } + inline const FlagGrid &getArg0() + { + return flags; + } + typedef FlagGrid type0; + void runMessage() + { + debMsg("Executing kernel CountEmptyCells ", 3); + debMsg("Kernel range" + << " x " << maxX << " y " << maxY << " z " << minZ << " - " << maxZ << " ", + 4); + }; + void operator()(const tbb::blocked_range<IndexInt> &__r) + { + for (IndexInt idx = __r.begin(); idx != (IndexInt)__r.end(); idx++) + op(idx, flags, numEmpty); + } + void run() + { + tbb::parallel_reduce(tbb::blocked_range<IndexInt>(0, size), *this); + } + CountEmptyCells(CountEmptyCells &o, tbb::split) : KernelBase(o), flags(o.flags), numEmpty(0) + { + } + void join(const CountEmptyCells &o) + { + numEmpty += o.numEmpty; + } + const FlagGrid &flags; + int numEmpty; +}; + +// ***************************************************************************** +// Misc helpers + +//! Change 'A' and 'rhs' such that pressure at 'fixPidx' is fixed to 'value' +void fixPressure(int fixPidx, + Real value, + Grid<Real> &rhs, + Grid<Real> &A0, + Grid<Real> &Ai, + Grid<Real> &Aj, + Grid<Real> &Ak) +{ + // Bring to rhs at neighbors + rhs[fixPidx + Ai.getStrideX()] -= Ai[fixPidx] * value; + rhs[fixPidx + Aj.getStrideY()] -= Aj[fixPidx] * value; + rhs[fixPidx - Ai.getStrideX()] -= Ai[fixPidx - Ai.getStrideX()] * value; + rhs[fixPidx - Aj.getStrideY()] -= Aj[fixPidx - Aj.getStrideY()] * value; + if (rhs.is3D()) { + rhs[fixPidx + Ak.getStrideZ()] -= Ak[fixPidx] * value; + rhs[fixPidx - Ak.getStrideZ()] -= Ak[fixPidx - Ak.getStrideZ()] * value; + } + + // Trivialize equation at 'fixPidx' to: pressure[fixPidx] = value + rhs[fixPidx] = value; + A0[fixPidx] = Real(1); + Ai[fixPidx] = Aj[fixPidx] = Ak[fixPidx] = Real(0); + Ai[fixPidx - Ai.getStrideX()] = Real(0); + Aj[fixPidx - Aj.getStrideY()] = Real(0); + if (rhs.is3D()) { + Ak[fixPidx - Ak.getStrideZ()] = Real(0); + } +} + +// for "static" MG mode, keep one MG data structure per fluid solver +// leave cleanup to OS/user if nonzero at program termination (PcMGStatic mode) +// alternatively, manually release in scene file with releaseMG +static std::map<FluidSolver *, GridMg *> gMapMG; + +void releaseMG(FluidSolver *solver = nullptr) +{ + // release all? + if (!solver) { + for (std::map<FluidSolver *, GridMg *>::iterator it = gMapMG.begin(); it != gMapMG.end(); + it++) { + if (it->first != nullptr) + releaseMG(it->first); + } + return; + } + + GridMg *mg = gMapMG[solver]; + if (mg) { + delete mg; + gMapMG[solver] = nullptr; + } +} +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, "releaseMG", !noTiming); + PyObject *_retval = 0; + { + ArgLocker _lock; + FluidSolver *solver = _args.getPtrOpt<FluidSolver>("solver", 0, nullptr, &_lock); + _retval = getPyNone(); + releaseMG(solver); + _args.check(); + } + pbFinalizePlugin(parent, "releaseMG", !noTiming); + return _retval; + } + catch (std::exception &e) { + pbSetError("releaseMG", e.what()); + return 0; + } +} +static const Pb::Register _RP_releaseMG("", "releaseMG", _W_0); +extern "C" { +void PbRegister_releaseMG() +{ + KEEP_UNUSED(_RP_releaseMG); +} +} + +// ***************************************************************************** +// Main pressure solve + +// Note , all three pressure solve helper functions take +// identical parameters, apart from the RHS grid (and different const values) + +//! Compute rhs for pressure solve + +void computePressureRhs(Grid<Real> &rhs, + const MACGrid &vel, + const Grid<Real> &pressure, + const FlagGrid &flags, + Real cgAccuracy = 1e-3, + const Grid<Real> *phi = 0, + const Grid<Real> *perCellCorr = 0, + const MACGrid *fractions = 0, + const MACGrid *obvel = 0, + Real gfClamp = 1e-04, + Real cgMaxIterFac = 1.5, + bool precondition = true, + int preconditioner = PcMIC, + bool enforceCompatibility = false, + bool useL2Norm = false, + bool zeroPressureFixing = false, + const Grid<Real> *curv = NULL, + const Real surfTens = 0.) +{ + // compute divergence and init right hand side + MakeRhs kernMakeRhs( + flags, rhs, vel, perCellCorr, fractions, obvel, phi, curv, surfTens, gfClamp); + + if (enforceCompatibility) + rhs += (Real)(-kernMakeRhs.sum / (Real)kernMakeRhs.cnt); +} +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, "computePressureRhs", !noTiming); + PyObject *_retval = 0; + { + ArgLocker _lock; + Grid<Real> &rhs = *_args.getPtr<Grid<Real>>("rhs", 0, &_lock); + const MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock); + const Grid<Real> &pressure = *_args.getPtr<Grid<Real>>("pressure", 2, &_lock); + const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 3, &_lock); + Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 4, 1e-3, &_lock); + const Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 5, 0, &_lock); + const Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>("perCellCorr", 6, 0, &_lock); + const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 7, 0, &_lock); + const MACGrid *obvel = _args.getPtrOpt<MACGrid>("obvel", 8, 0, &_lock); + Real gfClamp = _args.getOpt<Real>("gfClamp", 9, 1e-04, &_lock); + Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 10, 1.5, &_lock); + bool precondition = _args.getOpt<bool>("precondition", 11, true, &_lock); + int preconditioner = _args.getOpt<int>("preconditioner", 12, PcMIC, &_lock); + bool enforceCompatibility = _args.getOpt<bool>("enforceCompatibility", 13, false, &_lock); + bool useL2Norm = _args.getOpt<bool>("useL2Norm", 14, false, &_lock); + bool zeroPressureFixing = _args.getOpt<bool>("zeroPressureFixing", 15, false, &_lock); + const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 16, NULL, &_lock); + const Real surfTens = _args.getOpt<Real>("surfTens", 17, 0., &_lock); + _retval = getPyNone(); + computePressureRhs(rhs, + vel, + pressure, + flags, + cgAccuracy, + phi, + perCellCorr, + fractions, + obvel, + gfClamp, + cgMaxIterFac, + precondition, + preconditioner, + enforceCompatibility, + useL2Norm, + zeroPressureFixing, + curv, + surfTens); + _args.check(); + } + pbFinalizePlugin(parent, "computePressureRhs", !noTiming); + return _retval; + } + catch (std::exception &e) { + pbSetError("computePressureRhs", e.what()); + return 0; + } +} +static const Pb::Register _RP_computePressureRhs("", "computePressureRhs", _W_1); +extern "C" { +void PbRegister_computePressureRhs() +{ + KEEP_UNUSED(_RP_computePressureRhs); +} +} + +//! Build and solve pressure system of equations +//! perCellCorr: a divergence correction for each cell, optional +//! fractions: for 2nd order obstacle boundaries, optional +//! gfClamp: clamping threshold for ghost fluid method +//! cgMaxIterFac: heuristic to determine maximal number of CG iteations, increase for more accurate +//! solutions preconditioner: MIC, or MG (see Preconditioner enum) useL2Norm: use max norm by +//! default, can be turned to L2 here zeroPressureFixing: remove null space by fixing a single +//! pressure value, needed for MG curv: curvature for surface tension effects surfTens: surface +//! tension coefficient retRhs: return RHS divergence, e.g., for debugging; optional + +void solvePressureSystem(Grid<Real> &rhs, + MACGrid &vel, + Grid<Real> &pressure, + const FlagGrid &flags, + Real cgAccuracy = 1e-3, + const Grid<Real> *phi = 0, + const Grid<Real> *perCellCorr = 0, + const MACGrid *fractions = 0, + Real gfClamp = 1e-04, + Real cgMaxIterFac = 1.5, + bool precondition = true, + int preconditioner = PcMIC, + const bool enforceCompatibility = false, + const bool useL2Norm = false, + const bool zeroPressureFixing = false, + const Grid<Real> *curv = NULL, + const Real surfTens = 0.) +{ + if (precondition == false) + preconditioner = PcNone; // for backwards compatibility + + // reserve temp grids + FluidSolver *parent = flags.getParent(); + 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); + + // setup matrix and boundaries + MakeLaplaceMatrix(flags, A0, Ai, Aj, Ak, fractions); + + if (phi) { + ApplyGhostFluidDiagonal(A0, flags, *phi, gfClamp); + } + + // check whether we need to fix some pressure value... + // (manually enable, or automatically for high accuracy, can cause asymmetries otherwise) + if (zeroPressureFixing || cgAccuracy < 1e-07) { + if (FLOATINGPOINT_PRECISION == 1) + debMsg( + "Warning - high CG accuracy with single-precision floating point accuracy might not " + "converge...", + 2); + + int numEmpty = CountEmptyCells(flags); + IndexInt fixPidx = -1; + if (numEmpty == 0) { + // Determine appropriate fluid cell for pressure fixing + // 1) First check some preferred positions for approx. symmetric zeroPressureFixing + Vec3i topCenter( + flags.getSizeX() / 2, flags.getSizeY() - 1, flags.is3D() ? flags.getSizeZ() / 2 : 0); + Vec3i preferredPos[] = {topCenter, topCenter - Vec3i(0, 1, 0), topCenter - Vec3i(0, 2, 0)}; + + for (Vec3i pos : preferredPos) { + if (flags.isFluid(pos)) { + fixPidx = flags.index(pos); + break; + } + } + + // 2) Then search whole domain + if (fixPidx == -1) { + FOR_IJK_BND(flags, 1) + { + if (flags.isFluid(i, j, k)) { + fixPidx = flags.index(i, j, k); + // break FOR_IJK_BND loop + i = flags.getSizeX() - 1; + j = flags.getSizeY() - 1; + k = __kmax; + } + } + } + // debMsg("No empty cells! Fixing pressure of cell "<<fixPidx<<" to zero",1); + } + if (fixPidx >= 0) { + fixPressure(fixPidx, Real(0), rhs, A0, Ai, Aj, Ak); + static bool msgOnce = false; + if (!msgOnce) { + debMsg("Pinning pressure of cell " << fixPidx << " to zero", 2); + msgOnce = true; + } + } + } + + // CG setup + // note: the last factor increases the max iterations for 2d, which right now can't use a + // preconditioner + GridCgInterface *gcg; + if (vel.is3D()) + gcg = new GridCg<ApplyMatrix>(pressure, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak); + else + gcg = new GridCg<ApplyMatrix2D>( + pressure, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak); + + gcg->setAccuracy(cgAccuracy); + gcg->setUseL2Norm(useL2Norm); + + int maxIter = 0; + + Grid<Real> *pca0 = nullptr, *pca1 = nullptr, *pca2 = nullptr, *pca3 = nullptr; + GridMg *pmg = nullptr; + + // optional preconditioning + if (preconditioner == PcNone || preconditioner == PcMIC) { + maxIter = (int)(cgMaxIterFac * flags.getSize().max()) * (flags.is3D() ? 1 : 4); + + pca0 = new Grid<Real>(parent); + pca1 = new Grid<Real>(parent); + pca2 = new Grid<Real>(parent); + pca3 = new Grid<Real>(parent); + + gcg->setICPreconditioner(preconditioner == PcMIC ? GridCgInterface::PC_mICP : + GridCgInterface::PC_None, + pca0, + pca1, + pca2, + pca3); + } + else if (preconditioner == PcMGDynamic || preconditioner == PcMGStatic) { + maxIter = 100; + + pmg = gMapMG[parent]; + if (!pmg) { + pmg = new GridMg(pressure.getSize()); + gMapMG[parent] = pmg; + } + + gcg->setMGPreconditioner(GridCgInterface::PC_MGP, pmg); + } + + // CG solve + for (int iter = 0; iter < maxIter; iter++) { + if (!gcg->iterate()) + iter = maxIter; + if (iter < maxIter) + debMsg("FluidSolver::solvePressure iteration " << iter + << ", residual: " << gcg->getResNorm(), + 9); + } + debMsg("FluidSolver::solvePressure done. Iterations:" << gcg->getIterations() + << ", residual:" << gcg->getResNorm(), + 2); + + // Cleanup + if (gcg) + delete gcg; + if (pca0) + delete pca0; + if (pca1) + delete pca1; + if (pca2) + delete pca2; + if (pca3) + delete pca3; + + // PcMGDynamic: always delete multigrid solver after use + // PcMGStatic: keep multigrid solver for next solve + if (pmg && preconditioner == PcMGDynamic) + releaseMG(parent); +} +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, "solvePressureSystem", !noTiming); + PyObject *_retval = 0; + { + ArgLocker _lock; + Grid<Real> &rhs = *_args.getPtr<Grid<Real>>("rhs", 0, &_lock); + MACGrid &vel = *_args.getPtr<MACGrid>("vel", 1, &_lock); + Grid<Real> &pressure = *_args.getPtr<Grid<Real>>("pressure", 2, &_lock); + const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 3, &_lock); + Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 4, 1e-3, &_lock); + const Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 5, 0, &_lock); + const Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>("perCellCorr", 6, 0, &_lock); + const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 7, 0, &_lock); + Real gfClamp = _args.getOpt<Real>("gfClamp", 8, 1e-04, &_lock); + Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 9, 1.5, &_lock); + bool precondition = _args.getOpt<bool>("precondition", 10, true, &_lock); + int preconditioner = _args.getOpt<int>("preconditioner", 11, PcMIC, &_lock); + const bool enforceCompatibility = _args.getOpt<bool>( + "enforceCompatibility", 12, false, &_lock); + const bool useL2Norm = _args.getOpt<bool>("useL2Norm", 13, false, &_lock); + const bool zeroPressureFixing = _args.getOpt<bool>("zeroPressureFixing", 14, false, &_lock); + const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 15, NULL, &_lock); + const Real surfTens = _args.getOpt<Real>("surfTens", 16, 0., &_lock); + _retval = getPyNone(); + solvePressureSystem(rhs, + vel, + pressure, + flags, + cgAccuracy, + phi, + perCellCorr, + fractions, + gfClamp, + cgMaxIterFac, + precondition, + preconditioner, + enforceCompatibility, + useL2Norm, + zeroPressureFixing, + curv, + surfTens); + _args.check(); + } + pbFinalizePlugin(parent, "solvePressureSystem", !noTiming); + return _retval; + } + catch (std::exception &e) { + pbSetError("solvePressureSystem", e.what()); + return 0; + } +} +static const Pb::Register _RP_solvePressureSystem("", "solvePressureSystem", _W_2); +extern "C" { +void PbRegister_solvePressureSystem() +{ + KEEP_UNUSED(_RP_solvePressureSystem); +} +} + +//! Apply pressure gradient to make velocity field divergence free + +void correctVelocity(MACGrid &vel, + Grid<Real> &pressure, + const FlagGrid &flags, + Real cgAccuracy = 1e-3, + const Grid<Real> *phi = 0, + const Grid<Real> *perCellCorr = 0, + const MACGrid *fractions = 0, + Real gfClamp = 1e-04, + Real cgMaxIterFac = 1.5, + bool precondition = true, + int preconditioner = PcMIC, + bool enforceCompatibility = false, + bool useL2Norm = false, + bool zeroPressureFixing = false, + const Grid<Real> *curv = NULL, + const Real surfTens = 0.) +{ + knCorrectVelocity(flags, vel, pressure); + if (phi) { + knCorrectVelocityGhostFluid(vel, flags, pressure, *phi, gfClamp, curv, surfTens); + // improve behavior of clamping for large time steps: + knReplaceClampedGhostFluidVels(vel, flags, pressure, *phi, gfClamp); + } +} +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, "correctVelocity", !noTiming); + PyObject *_retval = 0; + { + ArgLocker _lock; + MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock); + Grid<Real> &pressure = *_args.getPtr<Grid<Real>>("pressure", 1, &_lock); + const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 2, &_lock); + Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 3, 1e-3, &_lock); + const Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 4, 0, &_lock); + const Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>("perCellCorr", 5, 0, &_lock); + const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 6, 0, &_lock); + Real gfClamp = _args.getOpt<Real>("gfClamp", 7, 1e-04, &_lock); + Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 8, 1.5, &_lock); + bool precondition = _args.getOpt<bool>("precondition", 9, true, &_lock); + int preconditioner = _args.getOpt<int>("preconditioner", 10, PcMIC, &_lock); + bool enforceCompatibility = _args.getOpt<bool>("enforceCompatibility", 11, false, &_lock); + bool useL2Norm = _args.getOpt<bool>("useL2Norm", 12, false, &_lock); + bool zeroPressureFixing = _args.getOpt<bool>("zeroPressureFixing", 13, false, &_lock); + const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 14, NULL, &_lock); + const Real surfTens = _args.getOpt<Real>("surfTens", 15, 0., &_lock); + _retval = getPyNone(); + correctVelocity(vel, + pressure, + flags, + cgAccuracy, + phi, + perCellCorr, + fractions, + gfClamp, + cgMaxIterFac, + precondition, + preconditioner, + enforceCompatibility, + useL2Norm, + zeroPressureFixing, + curv, + surfTens); + _args.check(); + } + pbFinalizePlugin(parent, "correctVelocity", !noTiming); + return _retval; + } + catch (std::exception &e) { + pbSetError("correctVelocity", e.what()); + return 0; + } +} +static const Pb::Register _RP_correctVelocity("", "correctVelocity", _W_3); +extern "C" { +void PbRegister_correctVelocity() +{ + KEEP_UNUSED(_RP_correctVelocity); +} +} + +//! Perform pressure projection of the velocity grid, calls +//! all three pressure helper functions in a row. + +void solvePressure(MACGrid &vel, + Grid<Real> &pressure, + const FlagGrid &flags, + Real cgAccuracy = 1e-3, + const Grid<Real> *phi = 0, + const Grid<Real> *perCellCorr = 0, + const MACGrid *fractions = 0, + const MACGrid *obvel = 0, + Real gfClamp = 1e-04, + Real cgMaxIterFac = 1.5, + bool precondition = true, + int preconditioner = PcMIC, + bool enforceCompatibility = false, + bool useL2Norm = false, + bool zeroPressureFixing = false, + const Grid<Real> *curv = NULL, + const Real surfTens = 0., + Grid<Real> *retRhs = NULL) +{ + Grid<Real> rhs(vel.getParent()); + + computePressureRhs(rhs, + vel, + pressure, + flags, + cgAccuracy, + phi, + perCellCorr, + fractions, + obvel, + gfClamp, + cgMaxIterFac, + precondition, + preconditioner, + enforceCompatibility, + useL2Norm, + zeroPressureFixing, + curv, + surfTens); + + solvePressureSystem(rhs, + vel, + pressure, + flags, + cgAccuracy, + phi, + perCellCorr, + fractions, + gfClamp, + cgMaxIterFac, + precondition, + preconditioner, + enforceCompatibility, + useL2Norm, + zeroPressureFixing, + curv, + surfTens); + + correctVelocity(vel, + pressure, + flags, + cgAccuracy, + phi, + perCellCorr, + fractions, + gfClamp, + cgMaxIterFac, + precondition, + preconditioner, + enforceCompatibility, + useL2Norm, + zeroPressureFixing, + curv, + surfTens); + + // optionally , return RHS + if (retRhs) { + retRhs->copyFrom(rhs); + } +} +static PyObject *_W_4(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, "solvePressure", !noTiming); + PyObject *_retval = 0; + { + ArgLocker _lock; + MACGrid &vel = *_args.getPtr<MACGrid>("vel", 0, &_lock); + Grid<Real> &pressure = *_args.getPtr<Grid<Real>>("pressure", 1, &_lock); + const FlagGrid &flags = *_args.getPtr<FlagGrid>("flags", 2, &_lock); + Real cgAccuracy = _args.getOpt<Real>("cgAccuracy", 3, 1e-3, &_lock); + const Grid<Real> *phi = _args.getPtrOpt<Grid<Real>>("phi", 4, 0, &_lock); + const Grid<Real> *perCellCorr = _args.getPtrOpt<Grid<Real>>("perCellCorr", 5, 0, &_lock); + const MACGrid *fractions = _args.getPtrOpt<MACGrid>("fractions", 6, 0, &_lock); + const MACGrid *obvel = _args.getPtrOpt<MACGrid>("obvel", 7, 0, &_lock); + Real gfClamp = _args.getOpt<Real>("gfClamp", 8, 1e-04, &_lock); + Real cgMaxIterFac = _args.getOpt<Real>("cgMaxIterFac", 9, 1.5, &_lock); + bool precondition = _args.getOpt<bool>("precondition", 10, true, &_lock); + int preconditioner = _args.getOpt<int>("preconditioner", 11, PcMIC, &_lock); + bool enforceCompatibility = _args.getOpt<bool>("enforceCompatibility", 12, false, &_lock); + bool useL2Norm = _args.getOpt<bool>("useL2Norm", 13, false, &_lock); + bool zeroPressureFixing = _args.getOpt<bool>("zeroPressureFixing", 14, false, &_lock); + const Grid<Real> *curv = _args.getPtrOpt<Grid<Real>>("curv", 15, NULL, &_lock); + const Real surfTens = _args.getOpt<Real>("surfTens", 16, 0., &_lock); + Grid<Real> *retRhs = _args.getPtrOpt<Grid<Real>>("retRhs", 17, NULL, &_lock); + _retval = getPyNone(); + solvePressure(vel, + pressure, + flags, + cgAccuracy, + phi, + perCellCorr, + fractions, + obvel, + gfClamp, + cgMaxIterFac, + precondition, + preconditioner, + enforceCompatibility, + useL2Norm, + zeroPressureFixing, + curv, + surfTens, + retRhs); + _args.check(); + } + pbFinalizePlugin(parent, "solvePressure", !noTiming); + return _retval; + } + catch (std::exception &e) { + pbSetError("solvePressure", e.what()); + return 0; + } +} +static const Pb::Register _RP_solvePressure("", "solvePressure", _W_4); +extern "C" { +void PbRegister_solvePressure() +{ + KEEP_UNUSED(_RP_solvePressure); +} +} + +} // namespace Manta |