Page MenuHome
Differential D3850 Diff 19899 extern/mantaflow/preprocessed/plugin/pressure.cpp

Changeset View

Standalone View

View Options

extern/mantaflow/preprocessed/plugin/pressure.cpp

  • This file was added.
// 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> &phi;
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> &phi;
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> &phi;
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