add vec grid diffusion for viscosity, also added LDC scene

scenes/lidDrivenCavity.py

@@ -0,0 +1,121 @@
+#
+# Lid driven cavity example with viscosity
+#
+from manta import *
+#setDebugLevel(2)
+
+# the normalized unit cube in manta has which world space size?
+worldScale = 1.0  
+
+# viscosity, in [m^2/s] , rescale to unit cube
+# uncomment one of these to select LDC with specific Reynolds nr
+# (higher ones will need larger resolution!)
+#visc       = 0.0002  / worldScale  # Re 5k
+visc       = 0.0001  / worldScale  # Re 10k
+#visc       = 0.00005 / worldScale  # Re 20k 
+#visc       = 0.00001 / worldScale  # Re 100k 
+#visc       = 0. # off, rely on numerical viscosity, no proper LDC!
+# move whole top side in one time unit, ie 1 m/s if domain is 1m 
+lidVel     = 1.00
+# reynolds nr , characteristic length = domain size 
+
+# to reduce the start up time, enable this to start the sim with a larger CFL number
+doQuickstart = True
+
+# compute Reynolds nr
+Re = 0.
+if visc>0.:
+	Re = lidVel * worldScale / visc
+
+# solver params
+res  = 100  
+gDim = 2
+gs = vec3(res,res,res)
+if gDim==2: gs.z = 1;
+s = Solver(name='main', gridSize = gs, dim=gDim)
+
+s.frameLength = 0.1 
+s.timestepMin = s.frameLength * 0.01
+s.timestepMax = s.frameLength * 1.0
+s.cfl         = 1.0
+s.timestep    = s.frameLength
+
+if doQuickstart:
+	# cfl reduction, start high, reduce over time
+	s.cfl         = 10.0 
+
+density  = s.create(RealGrid)
+# prepare grids
+diff1 = s.create(RealGrid)
+flags = s.create(FlagGrid)
+
+flags.initDomain(boundaryWidth=1) 
+flags.fillGrid()
+
+vel      = s.create(MACGrid)
+pressure = s.create(RealGrid)
+
+timings = Timings()
+
+if 1 and (GUI):
+	gui = Gui()
+	gui.show( True ) 
+	gui.pause()
+
+
+# geometry
+lid = s.create(Box, p0=gs*vec3(0.0,1.0-(1./float(gs.x)*3.1),0.0), p1=gs*vec3(1.0,1.0,1.0)  )
+# densities only for visual debugging... no influence!
+source = s.create(Cylinder, center=gs*vec3(0.5,0.5,0.5), radius=res*0.10, z=gs*vec3(0, 0.10, 0))
+
+
+# diffusion param for solve = const * dt / dx^2
+alphaV = visc * s.timestep * float(res*res)
+
+mantaMsg("Visc %f , a=%f , Re=%f " %(visc, alphaV, Re), 0 )
+
+
+#main loop
+lastFrame = -1
+outcnt    = 0
+for t in range(98765):
+	maxvel = vel.getMax()
+	s.adaptTimestep(maxvel)
+	mantaMsg('\nFrame %i, max vel %f, t %f, dt %f ' % (s.frame, maxvel, s.timeTotal, s.timestep ))
+
+	if doQuickstart:
+		# slowly reduce cfl number
+		if s.cfl>5.0 and s.timeTotal>20.:
+			s.cfl = 5.0
+		if s.cfl>1.0 and s.timeTotal>30.:
+			s.cfl = 1.0
+
+	lid.applyToGrid(grid=vel, value=Vec3( lidVel*float(gs.x),0,0) )
+	#if t%50==1:
+	if (lastFrame!=s.frame) and (s.frame%25==0):
+		source.applyToGrid(grid=density, value=1)
+
+	advectSemiLagrange(flags=flags, vel=vel, grid=density, order=2, clampMode=2) 
+	advectSemiLagrange(flags=flags, vel=vel, grid=vel,     order=2, clampMode=2)
+	resetOutflow(flags=flags,real=density) 
+
+	setWallBcs(flags=flags, vel=vel)    
+
+	# viscosity test 
+	density.setBound(0.0, 1)
+
+	# vel diffusion / viscosity!
+	if visc>0.:
+		cgSolveDiffusion( flags, vel, alphaV )
+
+
+	# solve pressure
+	setWallBcs(flags=flags, vel=vel)    
+	solvePressure(flags=flags, vel=vel, pressure=pressure)
+
+	if 0 and (GUI) and (lastFrame!=s.frame) and (s.frame%1==0):
+		gui.screenshot( 'ldc06re%05d_%04d.jpg' % (int(Re), s.frame) );
+
+	lastFrame = s.frame
+	s.step()
+

source/conjugategrad.cpp

@@ -17,7 +17,7 @@
 using namespace std;
 namespace Manta {
 
-const int CG_DEBUGLEVEL = 5;
+const int CG_DEBUGLEVEL = 2;
 
 //*****************************************************************************
 //  Precondition helpers
@@ -205,16 +205,18 @@ GridCg<APPLYMAT>::GridCg(Grid<Real>& dst, Grid<Real>& rhs, Grid<Real>& 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();            
+	//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) {
@@ -333,4 +335,92 @@ void GridCg<APPLYMAT>::setMGPreconditioner(PreconditionType method, GridMg* MG)
 template class GridCg<ApplyMatrix>;
 template class GridCg<ApplyMatrix2D>;
 
+
+
+//*****************************************************************************
+
+// diffusion , 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
+PYTHON() 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;
+}
+
+
+
+
 }; // DDF

source/conjugategrad.h

@@ -46,6 +46,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:
@@ -73,6 +76,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; }

source/grid.cpp

@@ -262,6 +262,31 @@ template<class T> std::string Grid<T>::getDataPointer() {
 
 // L1 / L2 functions
 
+/* NT_DEBUG merge 
+// L1
+KERNEL(idx, reduce=+) returns(double result=0.0)
+double knGridL1(const    Grid<Real>& a) { result += fabs(a[idx]); }
+
+KERNEL(idx, reduce=+) returns(double result=0.0)
+double knVecGridL1(const Grid<Vec3>& v) { result += norm(v[idx]); }
+
+//! compute L1 norm of whole grid content 
+PYTHON() Real getGridL1(Grid<Real>&    source) { return ( knGridL1   (source) ); }
+PYTHON() Real getVecGridL1(Grid<Vec3>& source) { return ( knVecGridL1(source) ); }
+
+// L2
+KERNEL(idx, reduce=+) returns(double result=0.0)
+double knGridL2(const    Grid<Real>& a) { result += a[idx]*a[idx]; }
+
+KERNEL(idx, reduce=+) returns(double result=0.0)
+double knVecGridL2(const Grid<Vec3>& v) { result += normSquare(v[idx]); }
+
+//! compute L2 norm of whole grid content 
+PYTHON() Real getGridL2(Grid<Real>&    source) { return sqrt( knGridL2   (source) ); }
+PYTHON() Real getVecGridL2(Grid<Vec3>& source) { return sqrt( knVecGridL2(source) ); }
+
+*/
+
 //! calculate L1 norm for whole grid with non-parallelized loop
 template<class GRID>
 Real loop_calcL1Grid (const GRID &grid, int bnd)
@@ -562,15 +587,15 @@ PYTHON() Real getGridAvg(Grid<Real>& source, FlagGrid* flags=NULL)
 
 //! transfer data between real and vec3 grids
 
-KERNEL(idx) void knGetComponent(Grid<Vec3>& source, Grid<Real>& target, int component) { 
+KERNEL(idx) void knGetComponent(const Grid<Vec3>& source, Grid<Real>& target, int component) { 
 	target[idx] = source[idx][component]; 
 }
-PYTHON() void getComponent(Grid<Vec3>& source, Grid<Real>& target, int component) { knGetComponent(source, target, component); }
+PYTHON() void getComponent(const Grid<Vec3>& source, Grid<Real>& target, int component) { knGetComponent(source, target, component); }
 
-KERNEL(idx) void knSetComponent(Grid<Real>& source, Grid<Vec3>& target, int component) { 
+KERNEL(idx) void knSetComponent(const Grid<Real>& source, Grid<Vec3>& target, int component) { 
 	target[idx][component] = source[idx]; 
 }
-PYTHON() void setComponent(Grid<Real>& source, Grid<Vec3>& target, int component) { knSetComponent(source, target, component); }
+PYTHON() void setComponent(const Grid<Real>& source, Grid<Vec3>& target, int component) { knSetComponent(source, target, component); }
 
 //******************************************************************************
 // Specialization classes

source/grid.h

@@ -351,8 +351,9 @@ void copyMacToVec3(MACGrid &source, Grid<Vec3>& target);
 void convertMacToVec3(MACGrid &source, Grid<Vec3>& target);
 void resampleVec3ToMac(Grid<Vec3>& source, MACGrid &target);
 void resampleMacToVec3 (MACGrid &source, Grid<Vec3>& target );
-void getComponent(const Grid<Vec3>& src, Grid<Real>& dst, int c);
-void setComponent(const Grid<Real>& src, Grid<Vec3>& dst, int c);
+
+void getComponent(const Grid<Vec3>& source, Grid<Real>& target, int component);
+void setComponent(const Grid<Real>& source, Grid<Vec3>& target, int component);