Page MenuHome
Differential D3850 Diff 17619 intern/mantaflow/intern/manta_develop/preprocessed/omp/plugin/pressure.cpp

Changeset View

Standalone View

View Options

intern/mantaflow/intern/manta_develop/preprocessed/omp/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 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),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 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 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 Grid<Real> * getArg5() {
return phi; }
typedef Grid<Real> type5;inline const Grid<Real> * getArg6() {
return curv; }
typedef Grid<Real> type6;inline const Real& getArg7() {
return surfTens; }
typedef Real type7;inline const Real& getArg8() {
return gfClamp; }
typedef Real type8; void runMessage() { debMsg("Executing kernel MakeRhs ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
int cnt = 0;double sum = 0;
#pragma omp for nowait
for (int k=minZ; k < maxZ; k++) for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,flags,rhs,vel,perCellCorr,fractions,phi,curv,surfTens,gfClamp,cnt,sum);
#pragma omp critical
{this->cnt += cnt; this->sum += sum; } }
}
else {
const int k=0;
#pragma omp parallel
{
int cnt = 0;double sum = 0;
#pragma omp for nowait
for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,flags,rhs,vel,perCellCorr,fractions,phi,curv,surfTens,gfClamp,cnt,sum);
#pragma omp critical
{this->cnt += cnt; this->sum += sum; } }
}
}
const FlagGrid& flags; Grid<Real>& rhs; const MACGrid& vel; const Grid<Real>* perCellCorr; const MACGrid* fractions; 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 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 run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
#pragma omp for
for (int k=minZ; k < maxZ; 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;
#pragma omp parallel
{
#pragma omp for
for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,flags,vel,pressure); }
}
}
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 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 run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
#pragma omp for
for (int k=minZ; k < maxZ; 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;
#pragma omp parallel
{
#pragma omp for
for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,A0,flags,phi,gfClamp); }
}
}
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 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 run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
#pragma omp for
for (int k=minZ; k < maxZ; 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;
#pragma omp parallel
{
#pragma omp for
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); }
}
}
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 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 run() {
const int _maxX = maxX; const int _maxY = maxY; if (maxZ > 1) {
#pragma omp parallel
{
#pragma omp for
for (int k=minZ; k < maxZ; 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;
#pragma omp parallel
{
#pragma omp for
for (int j=1; j < _maxY; j++) for (int i=1; i < _maxX; i++) op(i,j,k,vel,flags,pressure,phi,gfClamp); }
}
}
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 run() {
const IndexInt _sz = size;
#pragma omp parallel
{
int numEmpty = 0;
#pragma omp for nowait
for (IndexInt i = 0; i < _sz; i++) op(i,flags,numEmpty);
#pragma omp critical
{this->numEmpty += 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, 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, 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); 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); _retval = getPyNone(); computePressureRhs(rhs,vel,pressure,flags,cgAccuracy,phi,perCellCorr,fractions,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, 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, 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); 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); Grid<Real>* retRhs = _args.getPtrOpt<Grid<Real> >("retRhs",16,NULL,&_lock); _retval = getPyNone(); solvePressure(vel,pressure,flags,cgAccuracy,phi,perCellCorr,fractions,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); }
}
} // end namespace