Sebastián Barschkis
5539b68009
The solver will now automatically detect static scenes (no moving obstacles) and use a slightly faster pressure solve in those cases.
1517 lines
49 KiB
C++
1517 lines
49 KiB
C++
|
|
|
|
// 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];
|
|
// Release MG from previous step if present (e.g. if previous solve was with MGStatic)
|
|
if (pmg && preconditioner == PcMGDynamic) {
|
|
releaseMG(parent);
|
|
pmg = nullptr;
|
|
}
|
|
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
|