[Bf-blender-cvs] [dbc97a5cff3] fluid-mantaflow: updated manta pp files
Sebastián Barschkis
noreply at git.blender.org
Wed Nov 8 17:26:43 CET 2017
Commit: dbc97a5cff359cb446683d59f9af58084b371ff9
Author: Sebastián Barschkis
Date: Wed Nov 1 12:37:31 2017 +0100
Branches: fluid-mantaflow
https://developer.blender.org/rBdbc97a5cff359cb446683d59f9af58084b371ff9
updated manta pp files
===================================================================
M intern/mantaflow/CMakeLists.txt
M intern/mantaflow/intern/manta_pp/omp/conjugategrad.cpp
M intern/mantaflow/intern/manta_pp/omp/conjugategrad.h
M intern/mantaflow/intern/manta_pp/omp/fastmarch.cpp
M intern/mantaflow/intern/manta_pp/omp/fileio.cpp
M intern/mantaflow/intern/manta_pp/omp/fileio.h
M intern/mantaflow/intern/manta_pp/omp/fluidsolver.cpp
M intern/mantaflow/intern/manta_pp/omp/general.cpp
M intern/mantaflow/intern/manta_pp/omp/general.h
M intern/mantaflow/intern/manta_pp/omp/gitinfo.h
M intern/mantaflow/intern/manta_pp/omp/grid.cpp
M intern/mantaflow/intern/manta_pp/omp/grid.h
M intern/mantaflow/intern/manta_pp/omp/grid.h.reg
M intern/mantaflow/intern/manta_pp/omp/grid.h.reg.cpp
M intern/mantaflow/intern/manta_pp/omp/grid4d.h
M intern/mantaflow/intern/manta_pp/omp/grid4d.h.reg
M intern/mantaflow/intern/manta_pp/omp/grid4d.h.reg.cpp
M intern/mantaflow/intern/manta_pp/omp/noisefield.cpp
M intern/mantaflow/intern/manta_pp/omp/particle.cpp
M intern/mantaflow/intern/manta_pp/omp/particle.h
M intern/mantaflow/intern/manta_pp/omp/particle.h.reg
M intern/mantaflow/intern/manta_pp/omp/particle.h.reg.cpp
M intern/mantaflow/intern/manta_pp/omp/plugin/advection.cpp
A intern/mantaflow/intern/manta_pp/omp/plugin/apic.cpp
M intern/mantaflow/intern/manta_pp/omp/plugin/extforces.cpp
M intern/mantaflow/intern/manta_pp/omp/plugin/flip.cpp
M intern/mantaflow/intern/manta_pp/omp/plugin/fluidguiding.cpp
M intern/mantaflow/intern/manta_pp/omp/plugin/initplugins.cpp
M intern/mantaflow/intern/manta_pp/omp/plugin/pressure.cpp
A intern/mantaflow/intern/manta_pp/omp/pwrapper/numpyWrap.cpp
A intern/mantaflow/intern/manta_pp/omp/pwrapper/numpyWrap.h
M intern/mantaflow/intern/manta_pp/omp/pwrapper/pconvert.h
M intern/mantaflow/intern/manta_pp/omp/pwrapper/pvec3.cpp
M intern/mantaflow/intern/manta_pp/omp/pwrapper/pythonInclude.h
M intern/mantaflow/intern/manta_pp/omp/registration.cpp
M intern/mantaflow/intern/manta_pp/omp/test.cpp
M intern/mantaflow/intern/manta_pp/omp/turbulencepart.h.reg.cpp
M intern/mantaflow/intern/manta_pp/omp/util/randomstream.h
M intern/mantaflow/intern/manta_pp/omp/util/simpleimage.cpp
M intern/mantaflow/intern/manta_pp/omp/util/vector4d.cpp
M intern/mantaflow/intern/manta_pp/omp/util/vectorbase.cpp
M intern/mantaflow/intern/manta_pp/omp/util/vectorbase.h
M intern/mantaflow/intern/manta_pp/omp/vortexpart.h.reg.cpp
M intern/mantaflow/intern/manta_pp/tbb/conjugategrad.cpp
M intern/mantaflow/intern/manta_pp/tbb/conjugategrad.h
M intern/mantaflow/intern/manta_pp/tbb/fastmarch.cpp
M intern/mantaflow/intern/manta_pp/tbb/fileio.cpp
M intern/mantaflow/intern/manta_pp/tbb/fileio.h
M intern/mantaflow/intern/manta_pp/tbb/fluidsolver.cpp
M intern/mantaflow/intern/manta_pp/tbb/general.cpp
M intern/mantaflow/intern/manta_pp/tbb/general.h
M intern/mantaflow/intern/manta_pp/tbb/gitinfo.h
M intern/mantaflow/intern/manta_pp/tbb/grid.cpp
M intern/mantaflow/intern/manta_pp/tbb/grid.h
M intern/mantaflow/intern/manta_pp/tbb/grid.h.reg
M intern/mantaflow/intern/manta_pp/tbb/grid.h.reg.cpp
M intern/mantaflow/intern/manta_pp/tbb/grid4d.h
M intern/mantaflow/intern/manta_pp/tbb/grid4d.h.reg
M intern/mantaflow/intern/manta_pp/tbb/grid4d.h.reg.cpp
M intern/mantaflow/intern/manta_pp/tbb/noisefield.cpp
M intern/mantaflow/intern/manta_pp/tbb/particle.cpp
M intern/mantaflow/intern/manta_pp/tbb/particle.h
M intern/mantaflow/intern/manta_pp/tbb/particle.h.reg
M intern/mantaflow/intern/manta_pp/tbb/particle.h.reg.cpp
M intern/mantaflow/intern/manta_pp/tbb/plugin/advection.cpp
A intern/mantaflow/intern/manta_pp/tbb/plugin/apic.cpp
M intern/mantaflow/intern/manta_pp/tbb/plugin/extforces.cpp
M intern/mantaflow/intern/manta_pp/tbb/plugin/flip.cpp
M intern/mantaflow/intern/manta_pp/tbb/plugin/fluidguiding.cpp
M intern/mantaflow/intern/manta_pp/tbb/plugin/initplugins.cpp
M intern/mantaflow/intern/manta_pp/tbb/plugin/pressure.cpp
A intern/mantaflow/intern/manta_pp/tbb/pwrapper/numpyWrap.cpp
A intern/mantaflow/intern/manta_pp/tbb/pwrapper/numpyWrap.h
M intern/mantaflow/intern/manta_pp/tbb/pwrapper/pconvert.h
M intern/mantaflow/intern/manta_pp/tbb/pwrapper/pvec3.cpp
M intern/mantaflow/intern/manta_pp/tbb/pwrapper/pythonInclude.h
M intern/mantaflow/intern/manta_pp/tbb/registration.cpp
M intern/mantaflow/intern/manta_pp/tbb/test.cpp
M intern/mantaflow/intern/manta_pp/tbb/turbulencepart.h.reg.cpp
M intern/mantaflow/intern/manta_pp/tbb/util/randomstream.h
M intern/mantaflow/intern/manta_pp/tbb/util/simpleimage.cpp
M intern/mantaflow/intern/manta_pp/tbb/util/vector4d.cpp
M intern/mantaflow/intern/manta_pp/tbb/util/vectorbase.cpp
M intern/mantaflow/intern/manta_pp/tbb/util/vectorbase.h
M intern/mantaflow/intern/manta_pp/tbb/vortexpart.h.reg.cpp
===================================================================
diff --git a/intern/mantaflow/CMakeLists.txt b/intern/mantaflow/CMakeLists.txt
index e342a46364f..00df23dfa70 100644
--- a/intern/mantaflow/CMakeLists.txt
+++ b/intern/mantaflow/CMakeLists.txt
@@ -135,6 +135,7 @@ set(SRC
${MANTA_PP}/particle.h.reg
${MANTA_PP}/particle.h.reg.cpp
${MANTA_PP}/plugin/advection.cpp
+ ${MANTA_PP}/plugin/apic.cpp
${MANTA_PP}/plugin/extforces.cpp
${MANTA_PP}/plugin/fire.cpp
${MANTA_PP}/plugin/flip.cpp
@@ -142,12 +143,16 @@ set(SRC
${MANTA_PP}/plugin/initplugins.cpp
${MANTA_PP}/plugin/kepsilon.cpp
${MANTA_PP}/plugin/meshplugins.cpp
+# ${MANTA_PP}/plugin/numpyconvert.cpp
${MANTA_PP}/plugin/pressure.cpp
${MANTA_PP}/plugin/surfaceturbulence.cpp
+# ${MANTA_PP}/plugin/tfplugins.cpp
${MANTA_PP}/plugin/vortexplugins.cpp
${MANTA_PP}/plugin/waveletturbulence.cpp
${MANTA_PP}/plugin/waves.cpp
${MANTA_PP}/pwrapper/manta.h
+# ${MANTA_PP}/pwrapper/numpyWrap.cpp
+# ${MANTA_PP}/pwrapper/numpyWrap.h
${MANTA_PP}/pwrapper/pclass.cpp
${MANTA_PP}/pwrapper/pclass.h
${MANTA_PP}/pwrapper/pconvert.cpp
diff --git a/intern/mantaflow/intern/manta_pp/omp/conjugategrad.cpp b/intern/mantaflow/intern/manta_pp/omp/conjugategrad.cpp
index 68d5980fd78..f531391157c 100644
--- a/intern/mantaflow/intern/manta_pp/omp/conjugategrad.cpp
+++ b/intern/mantaflow/intern/manta_pp/omp/conjugategrad.cpp
@@ -19,7 +19,7 @@
* GNU General Public License (GPL)
* http://www.gnu.org/licenses
*
- * Conjugate gradient solver
+ * Conjugate gradient solver, for pressure and viscosity
*
******************************************************************************/
@@ -29,7 +29,7 @@
using namespace std;
namespace Manta {
-const int CG_DEBUGLEVEL = 5;
+const int CG_DEBUGLEVEL = 2;
//*****************************************************************************
// Precondition helpers
@@ -241,17 +241,14 @@ GridCg<APPLYMAT>::GridCg(Grid<Real>& dst, Grid<Real>& rhs, Grid<Real>& residual,
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)
-{
- dst.clear();
- residual.clear();
- search.clear();
- tmp.clear();
-}
+{ }
template<class APPLYMAT>
void GridCg<APPLYMAT>::doInit() {
mInited = true;
+ mIterations = 0;
+ mDst.clear();
mResidual.copyFrom( mRhs ); // p=0, residual = b
if (mPcMethod == PC_ICP) {
@@ -307,9 +304,8 @@ bool GridCg<APPLYMAT>::iterate() {
if(this->mUseL2Norm) {
mResNorm = GridSumSqr(mResidual).sum;
} else {
- mResNorm = mResidual.getMaxAbsValue();
+ mResNorm = mResidual.getMaxAbs();
}
- //if(mIterations % 10 == 9) debMsg("GridCg::Iteration i="<<mIterations<<", resNorm="<<mResNorm<<" accuracy="<<mAccuracy, 1);
// abort here to safe some work...
if(mResNorm<mAccuracy) {
@@ -325,6 +321,15 @@ bool GridCg<APPLYMAT>::iterate() {
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;
@@ -370,6 +375,93 @@ void GridCg<APPLYMAT>::setMGPreconditioner(PreconditionType method, GridMg* MG)
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(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; 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,&_loc [...]
+
+
+
+
}; // DDF
diff --git a/intern/mantaflow/intern/manta_pp/omp/conjugategrad.h b/intern/mantaflow/intern/manta_pp/omp/conjugategrad.h
index 1d7627cccd3..f93fd366c7b 100644
--- a/intern/mantaflow/intern/manta_pp/omp/conjugategrad.h
+++ b/intern/mantaflow/intern/manta_pp/omp/conjugategrad.h
@@ -58,6 +58,9 @@ class GridCgInterface {
virtual void setAccuracy(Real set) = 0;
virtual Real getAccuracy() const = 0;
+ //! force reinit upon next iterate() call, can be used for doing multiple solves
+ virtual void forceReinit() = 0;
+
void setUseL2Norm(bool set) { mUseL2Norm = set; }
protected:
@@ -85,6 +88,7 @@ class GridCg : public GridCgInterface {
//! init pointers, and copy values from "normal" matrix
void setICPreconditioner(PreconditionType method, Grid<Real> *A0, Grid<Real> *Ai, Grid<Real> *Aj, Grid<Real> *Ak);
void setMGPreconditioner(PreconditionType method, GridMg* MG);
+ void forceReinit() { mInited = false; }
// Accessors
Real getSigma() const { return mSigma; }
@@ -143,7 +147,7 @@ class GridCg : public GridCgInterface {
{
#pragma omp for
for (IndexInt i = 0; i < _sz; i++) op(i,flags,dst,src,A0,Ai,Aj,Ak); } } FlagGrid& flags; Grid<Real>& dst; Grid<Real>& src; Grid<Real>& A0; Grid<Real>& Ai; Grid<Real>& Aj; Grid<Real>& Ak; };
-#line 117 "conjugategrad.h"
+#line 121 "conjugategrad.h"
@@ -168,7 +172,7 @@ class GridCg : public GridCgInterface {
{
#pragma omp for
for (IndexInt i = 0; i < _sz; i++) op(i,flags,dst,src,A0,Ai,Aj,Ak); } } FlagGrid& flags; Grid<Real>& dst; Grid<Real>& src; Grid<Real>& A0; Grid<Real>& Ai; Grid<Real>& Aj; Grid<Real>& Ak; };
-#line 135 "conjugategrad.h"
+#line 139 "conjugategrad.h"
@@ -180,30 +184,30 @@ class GridCg : public GridCgInterface {
if(!fractions) {
// diagonal, A0
- if (!flags.isObstacle(i-1,j,k)) A0(i,j,k) += 1.;
- if (!flags.isObstacle(i+1,j,k)) A0(i,j,k) += 1.;
- if (!flags.isObstacle(i,j-1,k)) A0(i,j,k) += 1.;
- if (!flags.isObstacle(i,j+1,k)) A0(i,j,k) += 1.;
+ if (!flags.isObstacle(i-1,j,k))
@@ Diff output truncated at 10240 characters. @@
More information about the Bf-blender-cvs
mailing list