Adding vortex in cell, small optimization on examples

CHANGELOG.md

@@ -5,6 +5,15 @@ All notable changes to this project will be documented in this file.
 
 ### Added
 - Introduced getDomainIterator for Cell-list
+- Vortex in Cell example
+- Example to show how to add sensors in SPH/particle based methods (see)
+- Vortex in Cell example
+- Interpolation functions (see Numerics/vortex_in_cell example)
+- Gray-scott 3d example (see Grid/gray_scott_3d example)
+- HDF5 Check point restart for vector_dist particles (see ...) 
+- Raw reader for grid (see ...)
+- A way to specify names for proeprties and select properties to write
+- Ghost put on grid
 
 ### Fixed
 - Installation of PETSC in case with MUMPS try without MUMPS
@@ -12,6 +21,8 @@ All notable changes to this project will be documented in this file.
 - Bug in VTK writer binary in case of vectors
 - Bug in VTK writer binary: long int are not supported removing output
 
+
+
 ## [0.8.0] 28 February 2016
 
 ### Added

example/Grid/3_gray_scott_3d/Makefile

+include ../../example.mk
+
+CC=mpic++
+
+LDIR =
+
+OBJ = main.o
+
+%.o: %.cpp
+	$(CC) -O3 -c --std=c++11 -o $@ $< $(INCLUDE_PATH)
+
+gray_scott: $(OBJ)
+	$(CC) -o $@ $^ $(CFLAGS) $(LIBS_PATH) $(LIBS)
+
+all: gray_scott
+
+run: all
+	mpirun -np 4 ./gray_scott
+
+.PHONY: clean all run
+
+clean:
+	rm -f *.o *~ core gray_scott
+

example/Grid/3_gray_scott_3d/config.cfg

+[pack]
+files = main.cpp Makefile

example/Grid/3_gray_scott_3d/main.cpp

+#include "Grid/grid_dist_id.hpp"
+#include "data_type/aggregate.hpp"
+#include "timer.hpp"
+
+/*!
+ * \page Grid_3_gs Grid 3 Gray Scott in 3D
+ *
+ * # Solving a gray scott-system in 3D # {#e3_gs_gray_scott}
+ *
+ * This example is just an extension of the 2D Gray scott example.
+ * Here we show how to solve a non-linear reaction diffusion system in 3D
+ *
+ * \see \ref Grid_2_solve_eq
+ *
+ * \snippet Grid/3_gray_scott/main.cpp constants
+ * 
+ */
+
+//! \cond [constants] \endcond
+
+constexpr int U = 0;
+constexpr int V = 1;
+
+constexpr int x = 0;
+constexpr int y = 1;
+constexpr int z = 2;
+
+
+void init(grid_dist_id<3,double,aggregate<double,double> > & Old, grid_dist_id<3,double,aggregate<double,double> > & New, Box<3,double> & domain)
+{
+	auto it = Old.getDomainIterator();
+
+	while (it.isNext())
+	{
+		// Get the local grid key
+		auto key = it.get();
+
+		// Old values U and V
+		Old.template get<U>(key) = 1.0;
+		Old.template get<V>(key) = 0.0;
+
+		// Old values U and V
+		New.template get<U>(key) = 0.0;
+		New.template get<V>(key) = 0.0;
+
+		++it;
+	}
+
+	grid_key_dx<3> start({(long int)std::floor(Old.size(0)*1.55f/domain.getHigh(0)),(long int)std::floor(Old.size(1)*1.55f/domain.getHigh(1)),(long int)std::floor(Old.size(1)*1.55f/domain.getHigh(2))});
+	grid_key_dx<3> stop ({(long int)std::ceil (Old.size(0)*1.85f/domain.getHigh(0)),(long int)std::ceil (Old.size(1)*1.85f/domain.getHigh(1)),(long int)std::floor(Old.size(1)*1.85f/domain.getHigh(1))});
+	auto it_init = Old.getSubDomainIterator(start,stop);
+
+	while (it_init.isNext())
+	{
+		auto key = it_init.get();
+
+                Old.template get<U>(key) = 0.5 + (((double)std::rand())/RAND_MAX -0.5)/10.0;
+                Old.template get<V>(key) = 0.25 + (((double)std::rand())/RAND_MAX -0.5)/20.0;
+
+		++it_init;
+	}
+}
+
+//! \cond [end fun] \endcond
+
+
+int main(int argc, char* argv[])
+{
+	openfpm_init(&argc,&argv);
+
+	// domain
+	Box<3,double> domain({0.0,0.0},{2.5,2.5,2.5});
+	
+	// grid size
+	size_t sz[3] = {128,128,128};
+
+	// Define periodicity of the grid
+	periodicity<3> bc = {PERIODIC,PERIODIC,PERIODIC};
+	
+	// Ghost in grid unit
+	Ghost<3,long int> g(1);
+	
+	// deltaT
+	double deltaT = 1;
+
+	// Diffusion constant for specie U
+	double du = 2*1e-5;
+
+	// Diffusion constant for specie V
+	double dv = 1*1e-5;
+
+	// Number of timesteps
+        size_t timeSteps = 17000;
+
+	// K and F (Physical constant in the equation)
+        double K = 0.065;
+        double F = 0.034;
+
+	//! \cond [init lib] \endcond
+
+	/*!
+	 * \page Grid_3_gs Grid 3 Gray Scott
+	 *
+	 * Here we create 2 distributed grid in 2D Old and New. In particular because we want that
+	 * the second grid is distributed across processors in the same way we pass the decomposition
+	 * of the Old grid to the New one in the constructor with **Old.getDecomposition()**. Doing this,
+	 * we force the two grid to have the same decomposition.
+	 *
+	 * \snippet Grid/3_gray_scott/main.cpp init grid
+	 *
+	 */
+
+	//! \cond [init grid] \endcond
+
+	grid_dist_id<3, double, aggregate<double,double>> Old(sz,domain,g,bc);
+
+	// New grid with the decomposition of the old grid
+        grid_dist_id<3, double, aggregate<double,double>> New(Old.getDecomposition(),sz,g);
+
+	
+	// spacing of the grid on x and y
+	double spacing[3] = {Old.spacing(0),Old.spacing(1),Old.spacing(2)};
+
+	init(Old,New,domain);
+
+	// sync the ghost
+	size_t count = 0;
+	Old.template ghost_get<U,V>();
+
+	// because we assume that spacing[x] == spacing[y] we use formula 2
+	// and we calculate the prefactor of Eq 2
+	double uFactor = deltaT * du/(spacing[x]*spacing[x]);
+	double vFactor = deltaT * dv/(spacing[x]*spacing[x]);
+
+	for (size_t i = 0; i < timeSteps; ++i)
+	{
+		if (i % 300 == 0)
+			std::cout << "STEP: " << i << std::endl;
+
+		auto it = Old.getDomainIterator();
+
+		while (it.isNext())
+		{
+			auto key = it.get();
+
+			// update based on Eq 2
+			New.get<U>(key) = Old.get<U>(key) + uFactor * (
+										Old.get<U>(key.move(x,1)) +
+										Old.get<U>(key.move(x,-1)) +
+										Old.get<U>(key.move(y,1)) +
+										Old.get<U>(key.move(y,-1)) +
+										Old.get<U>(key.move(z,1)) +
+										Old.get<U>(key.move(z,-1)) -
+										6.0*Old.get<U>(key)) +
+										- deltaT * Old.get<U>(key) * Old.get<V>(key) * Old.get<V>(key) +
+										- deltaT * F * (Old.get<U>(key) - 1.0);
+
+
+			// update based on Eq 2
+			New.get<V>(key) = Old.get<V>(key) + vFactor * (
+										Old.get<V>(key.move(x,1)) +
+										Old.get<V>(key.move(x,-1)) +
+										Old.get<V>(key.move(y,1)) +
+										Old.get<V>(key.move(y,-1)) +
+										Old.get<V>(key.move(z,1)) +
+                                                                                Old.get<V>(key.move(z,-1)) -
+										6*Old.get<V>(key)) +
+										deltaT * Old.get<U>(key) * Old.get<V>(key) * Old.get<V>(key) +
+										- deltaT * (F+K) * Old.get<V>(key);
+
+			// Next point in the grid
+			++it;
+		}
+
+		// Here we copy New into the old grid in preparation of the new step
+		// It would be better to alternate, but using this we can show the usage
+		// of the function copy. To note that copy work only on two grid of the same
+		// decomposition. If you want to copy also the decomposition, or force to be
+		// exactly the same, use Old = New
+		Old.copy(New);
+
+		// After copy we synchronize again the ghost part U and V
+		Old.ghost_get<U,V>();
+
+		// Every 30 time step we output the configuration for
+		// visualization
+		if (i % 60 == 0)
+		{
+			Old.write("output",count,VTK_WRITER | FORMAT_BINARY);
+			count++;
+		}
+	}
+	
+	//! \cond [time stepping] \endcond
+
+	/*!
+	 * \page Grid_3_gs Grid 3 Gray Scott
+	 *
+	 * ## Finalize ##
+	 *
+	 * Deinitialize the library
+	 *
+	 * \snippet Grid/3_gray_scott/main.cpp finalize
+	 *
+	 */
+
+	//! \cond [finalize] \endcond
+
+	openfpm_finalize();
+
+	//! \cond [finalize] \endcond
+}

example/Numerics/Stoke_flow/1_3D_incompressible/main_petsc.cpp

@@ -465,7 +465,7 @@ int main(int argc, char* argv[])
 
 	//! \cond [copy write] \endcond
 
-	// Bring the solution to grid
+	// copy the solution to grid
 	fd.template copy<velocity,pressure>(x,{0,0},{sz[0]-1,sz[1]-1,sz[2]-1},g_dist);
 
 	g_dist.write("lid_driven_cavity_p_petsc");

example/Vector/7_SPH_dlb/main.cpp

@@ -324,30 +324,19 @@ inline void DWab(Point<3,double> & dx, Point<3,double> & DW, double r, bool prin
 {
 	const double qq=r/H;
 
-	if (qq < 1.0)
-	{
-		double qq2 = qq * qq;
-		double fac = (c1*qq + d1*qq2)/r;
+    double qq2 = qq * qq;
+    double fac1 = (c1*qq + d1*qq2)/r;
+    double b1 = (qq < 1.0)?1.0f:0.0f;
 
-		DW.get(0) = fac*dx.get(0);
-		DW.get(1) = fac*dx.get(1);
-		DW.get(2) = fac*dx.get(2);
-	}
-	else if (qq < 2.0)
-	{
-		double wqq = (2.0 - qq);
-		double fac = c2 * wqq * wqq / r;
+    double wqq = (2.0 - qq);
+    double fac2 = c2 * wqq * wqq / r;
+    double b2 = (qq >= 1.0 && qq < 2.0)?1.0f:0.0f;
 
-		DW.get(0) = fac * dx.get(0);
-		DW.get(1) = fac * dx.get(1);
-		DW.get(2) = fac * dx.get(2);
-	}
-	else
-	{
-		DW.get(0) = 0.0;
-		DW.get(1) = 0.0;
-		DW.get(2) = 0.0;
-	}
+    double factor = (b1*fac1 + b2*fac2);
+
+    DW.get(0) = factor * dx.get(0);
+    DW.get(1) = factor * dx.get(1);
+    DW.get(2) = factor * dx.get(2);
 }
 
 /*! \cond [kernel_sph_der] \endcond */

example/Vector/7_SPH_dlb_opt/main.cpp

@@ -209,10 +209,10 @@
  */
 
 // A constant to indicate boundary particles
-#define BOUNDARY 0
+#define BOUNDARY 1
 
 // A constant to indicate fluid particles
-#define FLUID 1
+#define FLUID 0
 
 // initial spacing between particles dp in the formulas
 const double dp = 0.0085;
@@ -229,6 +229,8 @@ const double gamma_ = 7.0;
 // sqrt(3.0*dp*dp) support of the kernel
 const double H = 0.0147224318643;
 
+const double FourH2 = 4.0 * H*H;
+
 // Eta in the formulas
 const double Eta2 = 0.01 * H*H;
 
@@ -280,47 +282,49 @@ double max_fluid_height = 0.0;
 // FLUID or BOUNDARY
 const size_t type = 0;
 
+// FLUID or BOUNDARY
+const size_t nn_num = 1;
+
 // Density
-const int rho = 1;
+const int rho = 2;
 
 // Density at step n-1
-const int rho_prev = 2;
+const int rho_prev = 3;
 
 // Pressure
-const int Pressure = 3;
+const int Pressure = 4;
 
 // Delta rho calculated in the force calculation
-const int drho = 4;
+const int drho = 5;
 
 // calculated force
-const int force = 5;
+const int force = 6;
 
 // velocity
-const int velocity = 6;
+const int velocity = 7;
 
 // velocity at previous step
-const int velocity_prev = 7;
+const int velocity_prev = 8;
 
 /*! \cond [sim parameters] \endcond */
 
 /*! \cond [vector_dist_def] \endcond */
 
 // Type of the vector containing particles
-typedef vector_dist<3,double,aggregate<size_t,double,  double,    double,     double,     double[3], double[3], double[3]> > particles;
+typedef vector_dist<3,double,aggregate<int, int,double,  double,    double,     double,     double[3], double[3], double[3]> > particles;
 //                                       |      |        |          |            |            |         |            |
 //                                       |      |        |          |            |            |         |            |
 //                                     type   density   density    Pressure    delta       force     velocity    velocity
 //                                                      at n-1                 density                           at n - 1
 
-
 struct ModelCustom
 {
 	template<typename Decomposition, typename vector> inline void addComputation(Decomposition & dec, vector & vd, size_t v, size_t p)
 	{
-		if (vd.template getProp<type>(p) == FLUID)
-			dec.addComputationCost(v,4);
-		else
-			dec.addComputationCost(v,3);
+                if (vd.template getProp<type>(p) == FLUID )
+                        dec.addComputationCost(v,4);
+                else
+                        dec.addComputationCost(v,3);
 	}
 
 	template<typename Decomposition> inline void applyModel(Decomposition & dec, size_t v)
@@ -328,13 +332,30 @@ struct ModelCustom
 		dec.setSubSubDomainComputationCost(v, dec.getSubSubDomainComputationCost(v) * dec.getSubSubDomainComputationCost(v));
 	}
 
-	double distributionTol()
+        float distributionTol()
 	{
 		return 1.01;
 	}
 };
 
 
+struct ModelCustom2
+{
+        template<typename Decomposition, typename vector> inline void addComputation(Decomposition & dec, vector & vd, size_t v, size_t p)
+        {
+            dec.addComputationCost(v,vd.template getProp<nn_num>(p) + 4);
+        }
+
+        template<typename Decomposition> inline void applyModel(Decomposition & dec, size_t v)
+        {
+        }
+
+        float distributionTol()
+        {
+                return 1.01;
+        }
+};
+
 inline void EqState(particles & vd)
 {
 	auto it = vd.getDomainIterator();
@@ -373,34 +394,23 @@ const double a2_4 = 0.25*a2;
 // Filled later
 double W_dap = 0.0;
 
-inline void DWab(Point<3,double> & dx, Point<3,double> & DW, double r, bool print)
+inline void DWab(Point<3,double> & dx, Point<3,double> & DW, double r)
 {
 	const double qq=r/H;
 
-	if (qq < 1.0)
-	{
-		double qq2 = qq * qq;
-		double fac = (c1*qq + d1*qq2)/r;
+    double qq2 = qq * qq;
+    double fac1 = (c1*qq + d1*qq2)/r;
+    double b1 = (qq < 1.0)?1.0f:0.0f;
 
-		DW.get(0) = fac*dx.get(0);
-		DW.get(1) = fac*dx.get(1);
-		DW.get(2) = fac*dx.get(2);
-	}
-	else if (qq < 2.0)
-	{
-		double wqq = (2.0 - qq);
-		double fac = c2 * wqq * wqq / r;
+    double wqq = (2.0 - qq);
+    double fac2 = c2 * wqq * wqq / r;
+    double b2 = (qq >= 1.0 && qq < 2.0)?1.0f:0.0f;
 
-		DW.get(0) = fac * dx.get(0);
-		DW.get(1) = fac * dx.get(1);
-		DW.get(2) = fac * dx.get(2);
-	}
-	else
-	{
-		DW.get(0) = 0.0;
-		DW.get(1) = 0.0;
-		DW.get(2) = 0.0;
-	}
+    double factor = (b1*fac1 + b2*fac2);
+
+    DW.get(0) = factor * dx.get(0);
+    DW.get(1) = factor * dx.get(1);
+    DW.get(2) = factor * dx.get(2);
 }
 
 inline double Tensile(double r, double rhoa, double rhob, double prs1, double prs2)
@@ -495,7 +505,7 @@ template<typename VerletList> inline double calc_forces(particles & vd, VerletLi
     /*! \cond [reset_particles2] \endcond */
 
     // Get an iterator over particles
-    auto part = vd.getParticleIteratorCRS(NN.getInternalCellList());
+    auto part = vd.getParticleIteratorCRS(NN);
 
 	double visc = 0;
 
@@ -508,6 +518,9 @@ template<typename VerletList> inline double calc_forces(particles & vd, VerletLi
 		// Get the position xp of the particle
 		Point<3,double> xa = vd.getPos(a);
 
+		// Type of the particle
+		size_t typea = vd.getProp<type>(a);
+
 		// Take the mass of the particle dependently if it is FLUID or BOUNDARY
 		double massa = (vd.getProp<type>(a) == FLUID)?MassFluid:MassBound;
 
@@ -520,149 +533,74 @@ template<typename VerletList> inline double calc_forces(particles & vd, VerletLi
 		// Get the Velocity of the particle a
 		Point<3,double> va = vd.getProp<velocity>(a);
 
-		// We threat FLUID particle differently from BOUNDARY PARTICLES ...
-		if (vd.getProp<type>(a) != FLUID)
-		{
-			// If it is a boundary particle calculate the delta rho based on equation 2
-			// This require to run across the neighborhoods particles of a
+        // Get an iterator over the neighborhood particles of p
+        auto Np = NN.template getNNIterator<NO_CHECK>(a);
 
-			//! \cond [nn_part] \endcond
-			auto Np = NN.template getNNIterator<NO_CHECK>(a);
-			//! \cond [nn_part] \endcond
+        size_t nn = 0;
 
-			// For each neighborhood particle
-			while (Np.isNext() == true)
-			{
-				// ... q
-				auto b = Np.get();
-
-				// Get the position xp of the particle
-				Point<3,double> xb = vd.getPos(b);
+        // For each neighborhood particle
+        while (Np.isNext() == true)
+        {
+                // ... q
+                auto b = Np.get();
 
-				//! \cond [skip_self] \endcond
+                // Get the position xp of the particle
+                Point<3,double> xb = vd.getPos(b);
 
-				// if (p == q) skip this particle
-				if (a == b)	{++Np; continue;};
+                // Get the distance between p and q
+                Point<3,double> dr = xa - xb;
+                // take the norm of this vector
+                float r2 = norm2(dr);
 
-				//! \cond [skip_self] \endcond
+                // if they interact
+                if (r2 < FourH2 && r2 > 1e-18)
+                {
+                        double r = sqrt(r2);
 
-				// get the mass of the particle
-				double massb = (vd.getProp<type>(b) == FLUID)?MassFluid:MassBound;
+                        unsigned int typeb = vd.getProp<type>(b);
 
-				// Get the velocity of the particle b
-				Point<3,double> vb = vd.getProp<velocity>(b);
+                        double massb = (typeb == FLUID)?MassFluid:MassBound;
+                        Point<3,double> vb = vd.getProp<velocity>(b);
+                        double Pb = vd.getProp<Pressure>(b);
+                        double rhob = vd.getProp<rho>(b);
 
-				// Get the pressure and density of particle b
-				double Pb = vd.getProp<Pressure>(b);
-				double rhob = vd.getProp<rho>(b);
+                        Point<3,double> v_rel = va - vb;
 
-				// Get the distance between p and q
-				Point<3,double> dr = xa - xb;
-				// take the norm of this vector
-				double r2 = norm2(dr);
+                        Point<3,double> DW;
+                        DWab(dr,DW,r);
 
-				// If the particles interact ...
-				if (r2 < 4.0*H*H)
-				{
-					// ... calculate delta rho
-					double r = sqrt(r2);
+                        //! \cond [symmetry2] \endcond
 
-					Point<3,double> dv = va - vb;
+                        double factor = - ((Pa + Pb) / (rhoa * rhob) + Tensile(r,rhoa,rhob,Pa,Pb) + Pi(dr,r2,v_rel,rhoa,rhob,massb,visc));
 
-					Point<3,double> DW;
-					DWab(dr,DW,r,false);
+                        // Bound - Bound does not produce any change
+                        factor = (typea == BOUNDARY && typeb == BOUNDARY)?0.0f:factor;
 
-					const double dot = dr.get(0)*dv.get(0) + dr.get(1)*dv.get(1) + dr.get(2)*dv.get(2);
-					const double dot_rr2 = dot/(r2+Eta2);
-					max_visc=std::max(dot_rr2,max_visc);
+                        vd.getProp<force>(a)[0] += massb * factor * DW.get(0);
+                        vd.getProp<force>(a)[1] += massb * factor * DW.get(1);
+                        vd.getProp<force>(a)[2] += massb * factor * DW.get(2);
 
-					double scal = (dv.get(0)*DW.get(0)+dv.get(1)*DW.get(1)+dv.get(2)*DW.get(2));