Big updates: add LeapFrog integration method (./Release/Sibernetic eu…

…ler or leapfrog), add posibility choose integratin method, add several posibilities define in comman arguments timestep of simulation and timelimit (duration of simulation in seconds) example ./Release/Sibernetic timestep=0.01 timelimit=1.

inc/owConfigProperty.h

@@ -37,12 +37,56 @@
 #define OWCONFIGURATION_H_
 
 #include "owOpenCLConstant.h"
+#include "owPhysicsConstant.h"
 
 struct owConfigProrerty{
 	//This value defines boundary of box in which simulation is
 	//Sizes of the box containing simulated 'world'
 	//Sizes choice is realized this way because it should be proportional to smoothing radius h
 public:
+	const int getParticleCount(){ return PARTICLE_COUNT; };
+	void setParticleCount(int value){
+		PARTICLE_COUNT = value;
+		PARTICLE_COUNT_RoundedUp = ((( PARTICLE_COUNT - 1 ) / local_NDRange_size ) + 1 ) * local_NDRange_size;
+	};
+	void setDeviceType(DEVICE type) { preferable_device_type = type; };
+	const int getParticleCount_RoundUp(){ return PARTICLE_COUNT_RoundedUp; };
+	const int getDeviceType() const { return preferable_device_type; };
+	const int getNumberOfIteration() const { return totalNumberOfIteration ;};
+	INTEGRATOR getIntegrationMethod() const { return integration_method; };
+	// Constructor
+	owConfigProrerty(int argc, char** argv){
+		preferable_device_type = CPU;
+		time_step = timeStep;
+		time_limit = 0.f;
+		beta = ::beta;
+		integration_method = EULER;
+		std::string s_temp;
+		for(int i = 1; i<argc; i++){
+			s_temp = argv[i];
+			if(s_temp.find("device=") == 0){
+				if(s_temp.find("GPU") != std::string::npos || s_temp.find("gpu") != std::string::npos)
+					preferable_device_type = GPU;
+			}
+			if(s_temp.find("timestep=") == 0){
+
+				time_step = ::atof( s_temp.substr(s_temp.find('=')+1).c_str());
+				time_step = (time_step > 0) ? time_step : timeStep;
+				//also we should recalculate beta if time_step is different from default value of timeStep in owPhysicsConstant
+				beta = time_step*time_step*mass*mass*2/(rho0*rho0);
+			}
+			if(s_temp.find("timelimit=") == 0){
+				time_limit = ::atof( s_temp.substr(s_temp.find('=')+1).c_str());
+			}
+			if(s_temp.find("LEAPFROG") != std::string::npos || s_temp.find("leapfrog") != std::string::npos){
+				integration_method = LEAPFROG;
+			}
+		}
+		totalNumberOfIteration = time_limit/time_step; // if it equals to 0 it means that simulation will work infinitely
+		calcDelta();
+	};
+	float getTimeStep() const { return timeStep; };
+	float getDelta() const { return delta; };
 	float xmin;
 	float xmax;
 	float ymin;
@@ -53,20 +97,68 @@ struct owConfigProrerty{
 	int gridCellsY;
 	int gridCellsZ;
 	int gridCellCount;
-	INTEGRATOR integration_method;
-	const int getParticleCount(){ return PARTICLE_COUNT; };
-	void setParticleCount(int value){
-		PARTICLE_COUNT = value;
-		PARTICLE_COUNT_RoundedUp = ((( PARTICLE_COUNT - 1 ) / local_NDRange_size ) + 1 ) * local_NDRange_size;
-	};
-	const int getParticleCount_RoundUp(){ return PARTICLE_COUNT_RoundedUp; };
-	void setDeviceType(DEVICE type){ preferable_device_type=type; };
-	const int getDeviceType(){ return preferable_device_type; };
 private:
+	/** Calculating delta parameter.
+	 *
+	 *  "In these situations,
+	 *  the SPH equations result in falsified values. To circumvent that problem, we pre-
+	 *  compute a single scaling factor δ according to the following formula [1, eq. 8] which is
+	 *  evaluated for a prototype particle with a filled neighborhood. The resulting value
+	 *  is then used for all particles. Finally, we end up with the following equations
+	 *  which are used in the PCISPH method" [1].
+	 *  [1] http://www.ifi.uzh.ch/vmml/publications/pcisph/pcisph.pdf
+	 */
+	inline void calcDelta(){
+
+		float x[] = { 1, 1, 0,-1,-1,-1, 0, 1, 1, 1, 0,-1,-1,-1, 0, 1, 1, 1, 0,-1,-1,-1, 0, 1, 2,-2, 0, 0, 0, 0, 0, 0 };
+	    float y[] = { 0, 1, 1, 1, 0,-1,-1,-1, 0, 1, 1, 1, 0,-1,-1,-1, 0, 1, 1, 1, 0,-1,-1,-1, 0, 0, 2,-2, 0, 0, 0, 0 };
+	    float z[] = { 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1,-1, 0, 0, 0, 0, 2,-2, 1,-1 };
+	    float sum1_x = 0.f;
+		float sum1_y = 0.f;
+		float sum1_z = 0.f;
+	    double sum1 = 0.0, sum2 = 0.0;
+		float v_x = 0.f;
+		float v_y = 0.f;
+		float v_z = 0.f;
+		float dist;
+		float particleRadius = pow(mass/rho0,1.f/3.f);  // It's equal to simulationScale
+														// TODO: replace it with simulation scale
+		float h_r_2;
+
+	    for (int i = 0; i < MAX_NEIGHBOR_COUNT; i++)
+	    {
+			v_x = x[i] * 0.8f/*1.f*/ * particleRadius; // return it back to 0.8 it's more stable
+			v_y = y[i] * 0.8f/*1.f*/ * particleRadius; // return it back to 0.8 it's more stable
+			v_z = z[i] * 0.8f/*1.f*/ * particleRadius; // return it back to 0.8 it's more stable
+
+	        dist = sqrt(v_x*v_x+v_y*v_y+v_z*v_z);//scaled, right?
+
+	        if (dist <= h*simulationScale)
+	        {
+				h_r_2 = pow((h*simulationScale - dist),2);//scaled
+
+	            sum1_x += h_r_2 * v_x / dist;
+				sum1_y += h_r_2 * v_y / dist;
+				sum1_z += h_r_2 * v_z / dist;
+
+	            sum2 += h_r_2 * h_r_2;
+	        }
+	    }
+		sum1 = sum1_x*sum1_x + sum1_y*sum1_y + sum1_z*sum1_z;
+		double result = 1.0 / (beta * gradWspikyCoefficient * gradWspikyCoefficient * (sum1 + sum2));
+		//return  1.0f / (beta * gradWspikyCoefficient * gradWspikyCoefficient * (sum1 + sum2));
+		delta = (float)result;
+	}
+
 	int PARTICLE_COUNT;
 	int PARTICLE_COUNT_RoundedUp;
+	int totalNumberOfIteration;
+	float time_step;
+	float time_limit;
+	float beta;
+	float delta;
 	DEVICE preferable_device_type;// 0-CPU, 1-GPU
+	INTEGRATOR integration_method; //DEFAULT is EULER
 };
 
-
 #endif /* OWCONFIGURATION_H_ */

inc/owPhysicTest.h

@@ -10,7 +10,7 @@
 
 #include "owPhysicsFluidSimulator.h"
 
-void test_energy_conservation();
+void test_energy_conservation(int argc, char **argv);
 
 
 #endif /* OWPHYSICTEST_H_ */

inc/owPhysicsConstant.h

@@ -130,7 +130,6 @@ const double divgradWviscosityCoefficient = - gradWspikyCoefficient;
 const float gravity_x = 0.0f;                       // Value of vector Gravity component x
 const float gravity_y = -9.8f;                      // Value of vector Gravity component y
 const float gravity_z = 0.0f;                       // Value of vector Gravity component z
-extern const float delta;                           // NOTE more info about this parameter see in file owPhysicsFluidSimulator.cpp description of function calcDelta()
 const int maxIteration = 3;                         // Number of iterations for Predictive-Corrective scheme
 
 const float mass_mult_Wpoly6Coefficient = (float) ( (double)mass * Wpoly6Coefficient );                       // Conversion of double value to float. For work with only 1st precision arithmetic.

inc/owPhysicsFluidSimulator.h

@@ -45,7 +45,7 @@ class owPhysicsFluidSimulator
 {
 public:
 	owPhysicsFluidSimulator(void);
-	owPhysicsFluidSimulator(owHelper * helper, DEVICE dev_type=CPU);
+	owPhysicsFluidSimulator(owHelper * helper, int argc, char ** argv);
 	~owPhysicsFluidSimulator(void);
 	/** Getter for position_cpp
 	 *
@@ -113,9 +113,9 @@ class owPhysicsFluidSimulator
 	void reset();
 private:
 	owOpenCLSolver * ocl_solver;
-	float * position_cpp;				// everywhere in the code %variableName%_cpp means that we create 
-	float * velocity_cpp;				// and initialize in 'ordinary' memory some data, which will be 
-	float * elasticConnectionsData_cpp; // copied later to OpenCL buffer %variableName% 
+	float * position_cpp;				// everywhere in the code %variableName%_cpp means that we create
+	float * velocity_cpp;				// and initialize in 'ordinary' memory some data, which will be
+	float * elasticConnectionsData_cpp; // copied later to OpenCL buffer %variableName%
 	int	  * membraneData_cpp;
 	int   * particleMembranesList_cpp;
 	//Helper arrays for displaying information about density changes

src/main.cpp

@@ -66,7 +66,7 @@ int main(int argc, char **argv)
 			}
 		}
 		if(run_tests){
-			test_energy_conservation();
+			test_energy_conservation(argc, argv);
 		}
 		else
 			run( argc, argv, graph, load_to );

src/owOpenCLSolver.cpp

@@ -226,14 +226,15 @@ void owOpenCLSolver::initializeOpenCL(owConfigProrerty * config)
     cl_device_type type;
 	const int device_type [] = {CL_DEVICE_TYPE_CPU,CL_DEVICE_TYPE_GPU};
 
-	unsigned int plList = 0;//selected platform index in platformList array [choose CPU by default]
+	int plList = -1;//selected platform index in platformList array [choose CPU by default]
 							//added autodetection of device number corresonding to preferrable device type (CPU|GPU) | otherwise the choice will be made from list of existing devices
 	cl_uint ciDeviceCount;
 	cl_device_id * devices_t;
 	bool bPassed = true, findDevice = false;
 	cl_int result;
 	cl_uint device_coumpute_unit_num;
 	cl_uint device_coumpute_unit_num_current = 0;
+	unsigned int deviceNum = 0;
 	//Selection of more appropriate device
 	for(int clSelectedPlatformID = 0;clSelectedPlatformID < (int)n_pl;clSelectedPlatformID++){
 		//if(findDevice)
@@ -253,14 +254,24 @@ void owOpenCLSolver::initializeOpenCL(owConfigProrerty * config)
 							plList = clSelectedPlatformID;
 							device_coumpute_unit_num_current = device_coumpute_unit_num;
 							findDevice = true;
+							deviceNum = i;
 						}
 						//break;
 					}
 				}
+				if(ciDeviceCount != 0 && plList < 0){
+					plList = clSelectedPlatformID;
+				}
 			}
 		}
 	}
-	if(!findDevice) plList = 0;
+	if(!findDevice){
+		//plList = 0;
+		deviceNum = 0;
+		std::cout << "Unfortunately OpenCL couldn't find device " << ((config->getDeviceType() == CPU)? "CPU":"GPU") << std::endl;
+		std::cout << "OpenCL try to init existing device " << std::endl;
+		config->setDeviceType((config->getDeviceType() == CPU)? GPU:CPU);
+	}
 	cl_context_properties cprops[3] = { CL_CONTEXT_PLATFORM, (cl_context_properties) (platformList[plList])(), 0 };
 	context = cl::Context( device_type[config->getDeviceType()], cprops, NULL, NULL, &err );
 	devices = context.getInfo< CL_CONTEXT_DEVICES >();
@@ -272,10 +283,11 @@ void owOpenCLSolver::initializeOpenCL(owConfigProrerty * config)
     unsigned long val2;
     size_t val3;
 	//uint deviceNum = 0;// causes "error C2065: 'uint' : undeclared identifier"
-    unsigned int deviceNum = 0;
-	result = devices[deviceNum].getInfo(CL_DEVICE_NAME,&cBuffer);// CL_INVALID_VALUE = -30;
+    result = devices[deviceNum].getInfo(CL_DEVICE_NAME,&cBuffer);// CL_INVALID_VALUE = -30;
 	if(result == CL_SUCCESS) std::cout << "CL_CONTEXT_PLATFORM ["<< plList << "]: CL_DEVICE_NAME [" << deviceNum << "]:\t" << cBuffer << "\n" << std::endl;
 	if(strlen(cBuffer)<1000) strcpy(device_full_name,cBuffer);
+	result = devices[deviceNum].getInfo(CL_DEVICE_TYPE,&cBuffer);
+	if(result == CL_SUCCESS) std::cout << "CL_CONTEXT_PLATFORM ["<< plList << "]: CL_DEVICE_TYPE [" << deviceNum << "]:\t" << (((int)cBuffer[0] == CL_DEVICE_TYPE_CPU)? "CPU" : "GPU") << std::endl;
 	result = devices[deviceNum].getInfo(CL_DEVICE_MAX_WORK_GROUP_SIZE,&val3);
 	if(result == CL_SUCCESS) std::cout << "CL_CONTEXT_PLATFORM ["<< plList << "]: CL_DEVICE_MAX_WORK_GROUP_SIZE [" <<  deviceNum <<"]: \t" << val3 <<std::endl;
 	result = devices[deviceNum].getInfo(CL_DEVICE_MAX_COMPUTE_UNITS,&value);
@@ -690,7 +702,7 @@ unsigned int owOpenCLSolver::_run_pcisph_predictPositions(owConfigProrerty * con
 	pcisph_predictPositions.setArg( 6, gravity_y );
 	pcisph_predictPositions.setArg( 7, gravity_z );
 	pcisph_predictPositions.setArg( 8, simulationScaleInv );
-	pcisph_predictPositions.setArg( 9, timeStep );
+	pcisph_predictPositions.setArg( 9, config->getTimeStep() );
 	pcisph_predictPositions.setArg(10, position );
 	pcisph_predictPositions.setArg(11, velocity );
 	pcisph_predictPositions.setArg(12, r0 );
@@ -760,7 +772,7 @@ unsigned int owOpenCLSolver::_run_pcisph_correctPressure(owConfigProrerty * conf
 	pcisph_correctPressure.setArg( 1, rho0 );
 	pcisph_correctPressure.setArg( 2, pressure );
 	pcisph_correctPressure.setArg( 3, rho );
-	pcisph_correctPressure.setArg( 4, delta );
+	pcisph_correctPressure.setArg( 4, config->getDelta() );
 	pcisph_correctPressure.setArg( 5, config->getParticleCount() );
 	int err = queue.enqueueNDRangeKernel(
 		pcisph_correctPressure, cl::NullRange, cl::NDRange( (int) (  config->getParticleCount_RoundUp() ) ),
@@ -793,7 +805,7 @@ unsigned int owOpenCLSolver::_run_pcisph_computePressureForceAcceleration(owConf
 	pcisph_computePressureForceAcceleration.setArg( 3, sortedPosition );
 	pcisph_computePressureForceAcceleration.setArg( 4, sortedVelocity );
 	pcisph_computePressureForceAcceleration.setArg( 5, particleIndexBack );
-	pcisph_computePressureForceAcceleration.setArg( 6, delta );
+	pcisph_computePressureForceAcceleration.setArg( 6, config->getDelta() );
 	pcisph_computePressureForceAcceleration.setArg( 7, mass_mult_gradWspikyCoefficient );
 	pcisph_computePressureForceAcceleration.setArg( 8, h );
 	pcisph_computePressureForceAcceleration.setArg( 9, simulationScale );
@@ -930,7 +942,7 @@ unsigned int owOpenCLSolver::_run_pcisph_integrate(int iterationCount, int pcisp
 	pcisph_integrate.setArg( 6, gravity_y );
 	pcisph_integrate.setArg( 7, gravity_z );
 	pcisph_integrate.setArg( 8, simulationScaleInv );
-	pcisph_integrate.setArg( 9, timeStep );
+	pcisph_integrate.setArg( 9, config->getTimeStep() );
 	pcisph_integrate.setArg( 10, config->xmin );
 	pcisph_integrate.setArg( 11, config->xmax );
 	pcisph_integrate.setArg( 12, config->ymin );

src/owPhysicsFluidSimulator.cpp

@@ -37,8 +37,6 @@
 
 #include "owPhysicsFluidSimulator.h"
 
-float calcDelta();
-extern const float delta = calcDelta();
 int numOfElasticConnections = 0; // Number of elastic connection TODO: move this to owConfig class
 int numOfLiquidP = 0;			 // Number of liquid particles TODO: move this to owConfig class
 int numOfElasticP = 0;			 // Number of liquid particles TODO: move this to owConfig class
@@ -55,15 +53,13 @@ int iter_step = 10;				 // Count of iteration which will be skipped before loggi
  *  @param dev_type
  *  defines preferable device type for current configuration
  */
-owPhysicsFluidSimulator::owPhysicsFluidSimulator(owHelper * helper,DEVICE dev_type)
+owPhysicsFluidSimulator::owPhysicsFluidSimulator(owHelper * helper,int argc, char ** argv)
 {
 	//int generateInitialConfiguration = 1;//1 to generate initial configuration, 0 - load from file
 
 	try{
 		iterationCount = 0;
-		config = new owConfigProrerty();
-		config->setDeviceType(dev_type);
-		config->integration_method = EULER;//LEAPFROG;
+		config = new owConfigProrerty(argc, argv);
 		// LOAD FROM FILE
 		owHelper::preLoadConfiguration(numOfMembranes, config, numOfLiquidP, numOfElasticP, numOfBoundaryP);
 
@@ -186,7 +182,7 @@ double owPhysicsFluidSimulator::simulationStep(const bool load_to)
 		ocl_solver->_runIndexPostPass(config);							helper->watch_report("_runIndexPostPass: \t%9.3f ms\n");
 		ocl_solver->_runFindNeighbors(config);							helper->watch_report("_runFindNeighbors: \t%9.3f ms\n");
 		//PCISPH PART
-		if(config->integration_method == LEAPFROG){ // in this case we should remmember value of position on stem i - 1
+		if(config->getIntegrationMethod() == LEAPFROG){ // in this case we should remmember value of position on stem i - 1
 			//Calc next time (t+dt) positions x(t+dt)
 			ocl_solver->_run_pcisph_integrate(iterationCount,0/*=positions_mode*/, config);
 		}
@@ -203,7 +199,7 @@ double owPhysicsFluidSimulator::simulationStep(const bool load_to)
 		}while( iter < maxIteration );
 
 		//and finally calculate v(t+dt)
-		if(config->integration_method == LEAPFROG){
+		if(config->getIntegrationMethod() == LEAPFROG){
 			ocl_solver->_run_pcisph_integrate(iterationCount,1/*=velocities_mode*/, config);		helper->watch_report("_runPCISPH: \t\t%9.3f ms\t3 iteration(s)\n");
 		}
 		else{
@@ -257,54 +253,3 @@ owPhysicsFluidSimulator::~owPhysicsFluidSimulator(void)
 	delete config;
 	delete ocl_solver;
 }
-/** Calculating delta parameter.
- *
- *  "In these situations,
- *  the SPH equations result in falsified values. To circumvent that problem, we pre-
- *  compute a single scaling factor δ according to the following formula [1, eq. 8] which is
- *  evaluated for a prototype particle with a filled neighborhood. The resulting value
- *  is then used for all particles. Finally, we end up with the following equations
- *  which are used in the PCISPH method" [1].
- *  [1] http://www.ifi.uzh.ch/vmml/publications/pcisph/pcisph.pdf
- */
-float calcDelta()
-{
-	float x[] = { 1, 1, 0,-1,-1,-1, 0, 1, 1, 1, 0,-1,-1,-1, 0, 1, 1, 1, 0,-1,-1,-1, 0, 1, 2,-2, 0, 0, 0, 0, 0, 0 };
-    float y[] = { 0, 1, 1, 1, 0,-1,-1,-1, 0, 1, 1, 1, 0,-1,-1,-1, 0, 1, 1, 1, 0,-1,-1,-1, 0, 0, 2,-2, 0, 0, 0, 0 };
-    float z[] = { 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1,-1,-1,-1,-1,-1,-1,-1,-1, 0, 0, 0, 0, 2,-2, 1,-1 };
-    float sum1_x = 0.f;
-	float sum1_y = 0.f;
-	float sum1_z = 0.f;
-    double sum1 = 0.0, sum2 = 0.0;
-	float v_x = 0.f;
-	float v_y = 0.f;
-	float v_z = 0.f;
-	float dist;
-	float particleRadius = pow(mass/rho0,1.f/3.f);  // It's equal to simulationScale
-													// TODO: replace it with simulation scale
-	float h_r_2;
-
-    for (int i = 0; i < MAX_NEIGHBOR_COUNT; i++)
-    {
-		v_x = x[i] * 0.8f/*1.f*/ * particleRadius; // return it back to 0.8 it's more stable
-		v_y = y[i] * 0.8f/*1.f*/ * particleRadius; // return it back to 0.8 it's more stable
-		v_z = z[i] * 0.8f/*1.f*/ * particleRadius; // return it back to 0.8 it's more stable
-
-        dist = sqrt(v_x*v_x+v_y*v_y+v_z*v_z);//scaled, right?
-
-        if (dist <= h*simulationScale)
-        {
-			h_r_2 = pow((h*simulationScale - dist),2);//scaled
-
-            sum1_x += h_r_2 * v_x / dist;
-			sum1_y += h_r_2 * v_y / dist;
-			sum1_z += h_r_2 * v_z / dist;
-
-            sum2 += h_r_2 * h_r_2;
-        }
-    }
-	sum1 = sum1_x*sum1_x + sum1_y*sum1_y + sum1_z*sum1_z;
-	double result = 1.0 / (beta * gradWspikyCoefficient * gradWspikyCoefficient * (sum1 + sum2));
-	//return  1.0f / (beta * gradWspikyCoefficient * gradWspikyCoefficient * (sum1 + sum2));
-	return (float)result;
-}