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. @@

_______________________________________________
Bf-blender-cvs mailing list
[email protected]
https://lists.blender.org/mailman/listinfo/bf-blender-cvs

Reply via email to