Page MenuHome
Differential D3850 Diff 17619 intern/mantaflow/intern/manta_develop/preprocessed/tbb/conjugategrad.cpp

Changeset View

Standalone View

View Options

intern/mantaflow/intern/manta_develop/preprocessed/tbb/conjugategrad.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
*
* Conjugate gradient solver, for pressure and viscosity
*
******************************************************************************/
#include "conjugategrad.h"
#include "commonkernels.h"
using namespace std;
namespace Manta {
const int CG_DEBUGLEVEL = 3;
//*****************************************************************************
// Precondition helpers
//! Preconditioning a la Wavelet Turbulence (needs 4 add. grids)
void InitPreconditionIncompCholesky(const FlagGrid& flags,
Grid<Real>& A0, Grid<Real>& Ai, Grid<Real>& Aj, Grid<Real>& Ak,
Grid<Real>& orgA0, Grid<Real>& orgAi, Grid<Real>& orgAj, Grid<Real>& orgAk)
{
// compute IC according to Golub and Van Loan
A0.copyFrom( orgA0 );
Ai.copyFrom( orgAi );
Aj.copyFrom( orgAj );
Ak.copyFrom( orgAk );
FOR_IJK(A0) {
if (flags.isFluid(i,j,k)) {
const IndexInt idx = A0.index(i,j,k);
A0[idx] = sqrt(A0[idx]);
// correct left and top stencil in other entries
// for i = k+1:n
// if (A(i,k) != 0)
// A(i,k) = A(i,k) / A(k,k)
Real invDiagonal = 1.0f / A0[idx];
Ai[idx] *= invDiagonal;
Aj[idx] *= invDiagonal;
Ak[idx] *= invDiagonal;
// correct the right and bottom stencil in other entries
// for j = k+1:n
// for i = j:n
// if (A(i,j) != 0)
// A(i,j) = A(i,j) - A(i,k) * A(j,k)
A0(i+1,j,k) -= square(Ai[idx]);
A0(i,j+1,k) -= square(Aj[idx]);
A0(i,j,k+1) -= square(Ak[idx]);
}
}
// invert A0 for faster computation later
InvertCheckFluid (flags, A0);
};
//! Preconditioning using modified IC ala Bridson (needs 1 add. grid)
void InitPreconditionModifiedIncompCholesky2(const FlagGrid& flags,
Grid<Real>&Aprecond,
Grid<Real>&A0, Grid<Real>& Ai, Grid<Real>& Aj, Grid<Real>& Ak)
{
// compute IC according to Golub and Van Loan
Aprecond.clear();
FOR_IJK(flags) {
if (!flags.isFluid(i,j,k)) continue;
const Real tau = 0.97;
const Real sigma = 0.25;
// compute modified incomplete cholesky
Real e = 0.;
e = A0(i,j,k)
- square(Ai(i-1,j,k) * Aprecond(i-1,j,k) )
- square(Aj(i,j-1,k) * Aprecond(i,j-1,k) )
- square(Ak(i,j,k-1) * Aprecond(i,j,k-1) ) ;
e -= tau * (
Ai(i-1,j,k) * ( Aj(i-1,j,k) + Ak(i-1,j,k) )* square( Aprecond(i-1,j,k) ) +
Aj(i,j-1,k) * ( Ai(i,j-1,k) + Ak(i,j-1,k) )* square( Aprecond(i,j-1,k) ) +
Ak(i,j,k-1) * ( Ai(i,j,k-1) + Aj(i,j,k-1) )* square( Aprecond(i,j,k-1) ) +
0. );
// stability cutoff
if(e < sigma * A0(i,j,k))
e = A0(i,j,k);
Aprecond(i,j,k) = 1. / sqrt( e );
}
};
//! Preconditioning using multigrid ala Dick et al.
void InitPreconditionMultigrid(GridMg* MG, Grid<Real>&A0, Grid<Real>& Ai, Grid<Real>& Aj, Grid<Real>& Ak, Real mAccuracy)
{
// build multigrid hierarchy if necessary
if (!MG->isASet()) MG->setA(&A0, &Ai, &Aj, &Ak);
MG->setCoarsestLevelAccuracy(mAccuracy * 1E-4);
MG->setSmoothing(1,1);
};
//! Apply WT-style ICP
void ApplyPreconditionIncompCholesky(Grid<Real>& dst, Grid<Real>& Var1, const FlagGrid& flags,
Grid<Real>& A0, Grid<Real>& Ai, Grid<Real>& Aj, Grid<Real>& Ak,
Grid<Real>& orgA0, Grid<Real>& orgAi, Grid<Real>& orgAj, Grid<Real>& orgAk)
{
// forward substitution
FOR_IJK(dst) {
if (!flags.isFluid(i,j,k)) continue;
dst(i,j,k) = A0(i,j,k) * (Var1(i,j,k)
- dst(i-1,j,k) * Ai(i-1,j,k)
- dst(i,j-1,k) * Aj(i,j-1,k)
- dst(i,j,k-1) * Ak(i,j,k-1));
}
// backward substitution
FOR_IJK_REVERSE(dst) {
const IndexInt idx = A0.index(i,j,k);
if (!flags.isFluid(idx)) continue;
dst[idx] = A0[idx] * ( dst[idx]
- dst(i+1,j,k) * Ai[idx]
- dst(i,j+1,k) * Aj[idx]
- dst(i,j,k+1) * Ak[idx]);
}
}
//! Apply Bridson-style mICP
void ApplyPreconditionModifiedIncompCholesky2(Grid<Real>& dst, Grid<Real>& Var1, const FlagGrid& flags,
Grid<Real>& Aprecond,
Grid<Real>& A0, Grid<Real>& Ai, Grid<Real>& Aj, Grid<Real>& Ak)
{
// forward substitution
FOR_IJK(dst) {
if (!flags.isFluid(i,j,k)) continue;
const Real p = Aprecond(i,j,k);
dst(i,j,k) = p * (Var1(i,j,k)
- dst(i-1,j,k) * Ai(i-1,j,k) * Aprecond(i-1,j,k)
- dst(i,j-1,k) * Aj(i,j-1,k) * Aprecond(i,j-1,k)
- dst(i,j,k-1) * Ak(i,j,k-1) * Aprecond(i,j,k-1) );
}
// backward substitution
FOR_IJK_REVERSE(dst) {
const IndexInt idx = A0.index(i,j,k);
if (!flags.isFluid(idx)) continue;
const Real p = Aprecond[idx];
dst[idx] = p * ( dst[idx]
- dst(i+1,j,k) * Ai[idx] * p
- dst(i,j+1,k) * Aj[idx] * p
- dst(i,j,k+1) * Ak[idx] * p);
}
}
//! Perform one Multigrid VCycle
void ApplyPreconditionMultigrid(GridMg* pMG, Grid<Real>& dst, Grid<Real>& Var1)
{
// one VCycle on "A*dst = Var1" with initial guess dst=0
pMG->setRhs(Var1);
pMG->doVCycle(dst);
}
//*****************************************************************************
// Kernels
//! Kernel: Compute the dot product between two Real grids
/*! Uses double precision internally */
struct GridDotProduct : public KernelBase {
GridDotProduct(const Grid<Real>& a, const Grid<Real>& b) : KernelBase(&a,0) ,a(a),b(b) ,result(0.0) {
runMessage(); run(); }
inline void op(IndexInt idx, const Grid<Real>& a, const Grid<Real>& b ,double& result) {
result += (a[idx] * b[idx]);
} inline operator double () {
return result; }
inline double & getRet() {
return result; }
inline const Grid<Real>& getArg0() {
return a; }
typedef Grid<Real> type0;inline const Grid<Real>& getArg1() {
return b; }
typedef Grid<Real> type1; void runMessage() { debMsg("Executing kernel GridDotProduct ", 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, a,b,result); }
void run() {
tbb::parallel_reduce (tbb::blocked_range<IndexInt>(0, size), *this); }
GridDotProduct (GridDotProduct& o, tbb::split) : KernelBase(o) ,a(o.a),b(o.b) ,result(0.0) {
}
void join(const GridDotProduct & o) {
result += o.result; }
const Grid<Real>& a; const Grid<Real>& b; double result; }
;;
//! Kernel: compute residual (init) and add to sigma
struct InitSigma : public KernelBase {
InitSigma(const FlagGrid& flags, Grid<Real>& dst, Grid<Real>& rhs, Grid<Real>& temp) : KernelBase(&flags,0) ,flags(flags),dst(dst),rhs(rhs),temp(temp) ,sigma(0) {
runMessage(); run(); }
inline void op(IndexInt idx, const FlagGrid& flags, Grid<Real>& dst, Grid<Real>& rhs, Grid<Real>& temp ,double& sigma) {
const double res = rhs[idx] - temp[idx];
dst[idx] = (Real)res;
// only compute residual in fluid region
if(flags.isFluid(idx))
sigma += res*res;
} inline operator double () {
return sigma; }
inline double & getRet() {
return sigma; }
inline const FlagGrid& getArg0() {
return flags; }
typedef FlagGrid type0;inline Grid<Real>& getArg1() {
return dst; }
typedef Grid<Real> type1;inline Grid<Real>& getArg2() {
return rhs; }
typedef Grid<Real> type2;inline Grid<Real>& getArg3() {
return temp; }
typedef Grid<Real> type3; void runMessage() { debMsg("Executing kernel InitSigma ", 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,dst,rhs,temp,sigma); }
void run() {
tbb::parallel_reduce (tbb::blocked_range<IndexInt>(0, size), *this); }
InitSigma (InitSigma& o, tbb::split) : KernelBase(o) ,flags(o.flags),dst(o.dst),rhs(o.rhs),temp(o.temp) ,sigma(0) {
}
void join(const InitSigma & o) {
sigma += o.sigma; }
const FlagGrid& flags; Grid<Real>& dst; Grid<Real>& rhs; Grid<Real>& temp; double sigma; }
;;
//! Kernel: update search vector
struct UpdateSearchVec : public KernelBase {
UpdateSearchVec(Grid<Real>& dst, Grid<Real>& src, Real factor) : KernelBase(&dst,0) ,dst(dst),src(src),factor(factor) {
runMessage(); run(); }
inline void op(IndexInt idx, Grid<Real>& dst, Grid<Real>& src, Real factor ) const {
dst[idx] = src[idx] + factor * dst[idx];
} inline Grid<Real>& getArg0() {
return dst; }
typedef Grid<Real> type0;inline Grid<Real>& getArg1() {
return src; }
typedef Grid<Real> type1;inline Real& getArg2() {
return factor; }
typedef Real type2; void runMessage() { debMsg("Executing kernel UpdateSearchVec ", 3); debMsg("Kernel range" << " x "<< maxX << " y "<< maxY << " z "<< minZ<<" - "<< maxZ << " " , 4); }; void operator() (const tbb::blocked_range<IndexInt>& __r) const {
for (IndexInt idx=__r.begin(); idx!=(IndexInt)__r.end(); idx++) op(idx, dst,src,factor); }
void run() {
tbb::parallel_for (tbb::blocked_range<IndexInt>(0, size), *this); }
Grid<Real>& dst; Grid<Real>& src; Real factor; }
;
//*****************************************************************************
// CG class
template<class APPLYMAT>
GridCg<APPLYMAT>::GridCg(Grid<Real>& dst, Grid<Real>& rhs, Grid<Real>& residual, Grid<Real>& search, const FlagGrid& flags, Grid<Real>& tmp,
Grid<Real>* pA0, Grid<Real>* pAi, Grid<Real>* pAj, Grid<Real>* pAk) :
GridCgInterface(), mInited(false), mIterations(0), mDst(dst), mRhs(rhs), mResidual(residual),
mSearch(search), mFlags(flags), mTmp(tmp), mpA0(pA0), mpAi(pAi), mpAj(pAj), mpAk(pAk),
mPcMethod(PC_None), mpPCA0(nullptr), mpPCAi(nullptr), mpPCAj(nullptr), mpPCAk(nullptr), mMG(nullptr), mSigma(0.), mAccuracy(VECTOR_EPSILON), mResNorm(1e20)
{ }
template<class APPLYMAT>
void GridCg<APPLYMAT>::doInit() {
mInited = true;
mIterations = 0;
mDst.clear();
mResidual.copyFrom( mRhs ); // p=0, residual = b
if (mPcMethod == PC_ICP) {
assertMsg(mDst.is3D(), "ICP only supports 3D grids so far");
InitPreconditionIncompCholesky(mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
ApplyPreconditionIncompCholesky(mTmp, mResidual, mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
} else if (mPcMethod == PC_mICP) {
assertMsg(mDst.is3D(), "mICP only supports 3D grids so far");
InitPreconditionModifiedIncompCholesky2(mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk);
ApplyPreconditionModifiedIncompCholesky2(mTmp, mResidual, mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk);
} else if (mPcMethod == PC_MGP) {
InitPreconditionMultigrid(mMG, *mpA0, *mpAi, *mpAj, *mpAk, mAccuracy);
ApplyPreconditionMultigrid(mMG, mTmp, mResidual);
} else {
mTmp.copyFrom( mResidual );
}
mSearch.copyFrom( mTmp );
mSigma = GridDotProduct(mTmp, mResidual);
}
template<class APPLYMAT>
bool GridCg<APPLYMAT>::iterate() {
if(!mInited) doInit();
mIterations++;
// create matrix application operator passed as template argument,
// this could reinterpret the mpA pointers (not so clean right now)
// tmp = applyMat(search)
APPLYMAT (mFlags, mTmp, mSearch, *mpA0, *mpAi, *mpAj, *mpAk);
// alpha = sigma/dot(tmp, search)
Real dp = GridDotProduct(mTmp, mSearch);
Real alpha = 0.;
if(fabs(dp)>0.) alpha = mSigma / (Real)dp;
gridScaledAdd<Real,Real>(mDst, mSearch, alpha); // dst += search * alpha
gridScaledAdd<Real,Real>(mResidual, mTmp, -alpha); // residual += tmp * -alpha
if (mPcMethod == PC_ICP)
ApplyPreconditionIncompCholesky(mTmp, mResidual, mFlags, *mpPCA0, *mpPCAi, *mpPCAj, *mpPCAk, *mpA0, *mpAi, *mpAj, *mpAk);
else if (mPcMethod == PC_mICP)
ApplyPreconditionModifiedIncompCholesky2(mTmp, mResidual, mFlags, *mpPCA0, *mpA0, *mpAi, *mpAj, *mpAk);
else if (mPcMethod == PC_MGP)
ApplyPreconditionMultigrid(mMG, mTmp, mResidual);
else
mTmp.copyFrom( mResidual );
// use the l2 norm of the residual for convergence check? (usually max norm is recommended instead)
if(this->mUseL2Norm) {
mResNorm = GridSumSqr(mResidual).sum;
} else {
mResNorm = mResidual.getMaxAbs();
}
// abort here to safe some work...
if(mResNorm<mAccuracy) {
mSigma = mResNorm; // this will be returned later on to the caller...
return false;
}
Real sigmaNew = GridDotProduct(mTmp, mResidual);
Real beta = sigmaNew / mSigma;
// search = tmp + beta * search
UpdateSearchVec (mSearch, mTmp, beta);
debMsg("GridCg::iterate i="<<mIterations<<" sigmaNew="<<sigmaNew<<" sigmaLast="<<mSigma<<" alpha="<<alpha<<" beta="<<beta<<" ", CG_DEBUGLEVEL);
mSigma = sigmaNew;
if(!(mResNorm<1e35)) {
if(mPcMethod == PC_MGP) {
// diverging solves can be caused by the static multigrid mode, we cannot detect this here, though
// only the pressure solve call "knows" whether the MG is static or dynamics...
debMsg("GridCg::iterate: Warning - this diverging solve can be caused by the 'static' mode of the MG preconditioner. If the static mode is active, try switching to dynamic.", 1);
}
errMsg("GridCg::iterate: The CG solver diverged, residual norm > 1e30, stopping.");
}
//debMsg("PB-CG-Norms::p"<<sqrt( GridOpNormNosqrt(mpDst, mpFlags).getValue() ) <<" search"<<sqrt( GridOpNormNosqrt(mpSearch, mpFlags).getValue(), CG_DEBUGLEVEL ) <<" res"<<sqrt( GridOpNormNosqrt(mpResidual, mpFlags).getValue() ) <<" tmp"<<sqrt( GridOpNormNosqrt(mpTmp, mpFlags).getValue() ), CG_DEBUGLEVEL ); // debug
return true;
}
template<class APPLYMAT>
void GridCg<APPLYMAT>::solve(int maxIter) {
for (int iter=0; iter<maxIter; iter++) {
if (!iterate()) iter=maxIter;
}
return;
}
static bool gPrint2dWarning = true;
template<class APPLYMAT>
void GridCg<APPLYMAT>::setICPreconditioner(PreconditionType method, Grid<Real> *A0, Grid<Real> *Ai, Grid<Real> *Aj, Grid<Real> *Ak) {
assertMsg(method==PC_ICP || method==PC_mICP, "GridCg<APPLYMAT>::setICPreconditioner: Invalid method specified.");
mPcMethod = method;
if( (!A0->is3D())) {
if(gPrint2dWarning) {
debMsg("ICP/mICP pre-conditioning only supported in 3D for now, disabling it.", 1);
gPrint2dWarning = false;
}
mPcMethod=PC_None;
}
mpPCA0 = A0;
mpPCAi = Ai;
mpPCAj = Aj;
mpPCAk = Ak;
}
template<class APPLYMAT>
void GridCg<APPLYMAT>::setMGPreconditioner(PreconditionType method, GridMg* MG) {
assertMsg(method==PC_MGP, "GridCg<APPLYMAT>::setMGPreconditioner: Invalid method specified.");
mPcMethod = method;
mMG = MG;
}
// explicit instantiation
template class GridCg<ApplyMatrix>;
template class GridCg<ApplyMatrix2D>;
//*****************************************************************************
// diffusion for real and vec grids, e.g. for viscosity
//! do a CG solve for diffusion; note: diffusion coefficient alpha given in grid space,
// rescale in python file for discretization independence (or physical correspondence)
// see lidDrivenCavity.py for an example
void cgSolveDiffusion(const FlagGrid& flags, GridBase& grid, Real alpha = 0.25, Real cgMaxIterFac = 1.0, Real cgAccuracy = 1e-4 ) {
// reserve temp grids
FluidSolver* parent = flags.getParent();
Grid<Real> rhs(parent);
Grid<Real> residual(parent), search(parent), tmp(parent);
Grid<Real> A0(parent), Ai(parent), Aj(parent), Ak(parent);
// setup matrix and boundaries
FlagGrid flagsDummy(parent);
flagsDummy.setConst(FlagGrid::TypeFluid);
MakeLaplaceMatrix (flagsDummy, A0, Ai, Aj, Ak);
FOR_IJK(flags) {
if(flags.isObstacle(i,j,k)) {
Ai(i,j,k) = Aj(i,j,k) = Ak(i,j,k) = 0.0;
A0(i,j,k) = 1.0;
} else {
Ai(i,j,k) *= alpha;
Aj(i,j,k) *= alpha;
Ak(i,j,k) *= alpha;
A0(i,j,k) *= alpha;
A0(i,j,k) += 1.;
}
}
GridCgInterface *gcg;
// note , no preconditioning for now...
const int maxIter = (int)(cgMaxIterFac * flags.getSize().max()) * (flags.is3D() ? 1 : 4);
if (grid.getType() & GridBase::TypeReal) {
Grid<Real>& u = ((Grid<Real>&) grid);
rhs.copyFrom(u);
if (flags.is3D())
gcg = new GridCg<ApplyMatrix >(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak );
else
gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak );
gcg->setAccuracy( cgAccuracy );
gcg->solve(maxIter);
debMsg("FluidSolver::solveDiffusion iterations:"<<gcg->getIterations()<<", res:"<<gcg->getSigma(), CG_DEBUGLEVEL);
}
else
if( (grid.getType() & GridBase::TypeVec3) || (grid.getType() & GridBase::TypeVec3) )
{
Grid<Vec3>& vec = ((Grid<Vec3>&) grid);
Grid<Real> u(parent);
// core solve is same as for a regular real grid
if (flags.is3D())
gcg = new GridCg<ApplyMatrix >(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak );
else
gcg = new GridCg<ApplyMatrix2D>(u, rhs, residual, search, flags, tmp, &A0, &Ai, &Aj, &Ak );
gcg->setAccuracy( cgAccuracy );
// diffuse every component separately
for(int component = 0; component< (grid.is3D() ? 3:2); ++component) {
getComponent( vec, u, component );
gcg->forceReinit();
rhs.copyFrom(u);
gcg->solve(maxIter);
debMsg("FluidSolver::solveDiffusion vec3, iterations:"<<gcg->getIterations()<<", res:"<<gcg->getSigma(), CG_DEBUGLEVEL);
setComponent( u, vec, component );
}
} else {
errMsg("cgSolveDiffusion: Grid Type is not supported (only Real, Vec3, MAC, or Levelset)");
}
delete gcg;
} 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, "cgSolveDiffusion" , !noTiming ); PyObject *_retval = 0; {
ArgLocker _lock; const FlagGrid& flags = *_args.getPtr<FlagGrid >("flags",0,&_lock); GridBase& grid = *_args.getPtr<GridBase >("grid",1,&_lock); Real alpha = _args.getOpt<Real >("alpha",2,0.25,&_lock); Real cgMaxIterFac = _args.getOpt<Real >("cgMaxIterFac",3,1.0,&_lock); Real cgAccuracy = _args.getOpt<Real >("cgAccuracy",4,1e-4 ,&_lock); _retval = getPyNone(); cgSolveDiffusion(flags,grid,alpha,cgMaxIterFac,cgAccuracy); _args.check(); }
pbFinalizePlugin(parent,"cgSolveDiffusion", !noTiming ); return _retval; }
catch(std::exception& e) {
pbSetError("cgSolveDiffusion",e.what()); return 0; }
}
static const Pb::Register _RP_cgSolveDiffusion ("","cgSolveDiffusion",_W_0);
extern "C" {
void PbRegister_cgSolveDiffusion() {
KEEP_UNUSED(_RP_cgSolveDiffusion); }
}
}; // DDF