Add LeapFrog Integration

inc/owConfigProperty.h

@@ -53,18 +53,19 @@ 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(int type){ preferable_device_type=type; };
+	void setDeviceType(DEVICE type){ preferable_device_type=type; };
 	const int getDeviceType(){ return preferable_device_type; };
 private:
 	int PARTICLE_COUNT;
 	int PARTICLE_COUNT_RoundedUp;
-	int preferable_device_type;// 0-CPU, 1-GPU
+	DEVICE preferable_device_type;// 0-CPU, 1-GPU
 };
 
 

inc/owOpenCLConstant.h

@@ -53,5 +53,6 @@
 const int local_NDRange_size = 256;
 
 enum DEVICE { CPU = 0, GPU = 1};
+enum INTEGRATOR { EULER = 0, LEAPFROG = 1 };
 
 #endif // #ifndef OW_OPENCL_CONSTANT_H

inc/owOpenCLSolver.h

@@ -88,14 +88,14 @@ class owOpenCLSolver
 	unsigned int _run_pcisph_predictDensity(owConfigProrerty * config);
 	unsigned int _run_pcisph_correctPressure(owConfigProrerty * config);
 	unsigned int _run_pcisph_computePressureForceAcceleration(owConfigProrerty * config);
-	unsigned int _run_pcisph_integrate(int iterationCount, owConfigProrerty * config);
+	unsigned int _run_pcisph_integrate(int iterationCount, int pcisph_integrate_mode, owConfigProrerty * config);
 	//Kernels for membrane handling interaction
 	unsigned int _run_clearMembraneBuffers(owConfigProrerty * config);
 	unsigned int _run_computeInteractionWithMembranes(owConfigProrerty * config);
 	unsigned int _run_computeInteractionWithMembranes_finalize(owConfigProrerty * config);
 	//
 	unsigned int updateMuscleActivityData(float *_muscle_activation_signal_cpp);
-	
+
 	void read_position_buffer( float * position_cpp, owConfigProrerty * config) { copy_buffer_from_device( position_cpp, position, config->getParticleCount() * sizeof( float ) * 4 ); };
 	void read_velocity_buffer( float * velocity_cpp, owConfigProrerty * config) { copy_buffer_from_device( velocity_cpp, velocity, config->getParticleCount() * sizeof( float ) * 4 ); };
 	void read_density_buffer( float * density_cpp, owConfigProrerty * config ) { copy_buffer_from_device( density_cpp, rho, config->getParticleCount() * sizeof( float ) * 1 ); }; // This need only for visualization current density of particle (graphic effect)
@@ -111,7 +111,7 @@ class owOpenCLSolver
 	cl::CommandQueue		  queue;
 	cl::Program				  program;
 	// Buffers
-	cl::Buffer  muscle_activation_signal;   // array storing data (activation signals) for an array of muscles. 
+	cl::Buffer  muscle_activation_signal;   // array storing data (activation signals) for an array of muscles.
                                             // now each can be activated by user independently
 
 	cl::Buffer acceleration;                // Acceleration buffer

src/owHelper.cpp

@@ -87,7 +87,8 @@ void owHelper::refreshTime()
 #endif
 }
 
-//READ DEFAULT CONFIGURATATION FROM FILE IN CONFIGURATION FOLDER
+
+//READ DEFAULT CONFIGURATATION FROM FILE IN CONFIGURATION FOLDER TODO move it into configuration struct
 int read_position = 0;
 std::string owHelper::path = "./configuration/";
 std::string owHelper::suffix = "";
@@ -526,11 +527,11 @@ void owHelper::watch_report( const char * str )
 #elif defined(__APPLE__)
     uint64_t elapsedNano;
     static mach_timebase_info_data_t    sTimebaseInfo;
-    
+
     if ( sTimebaseInfo.denom == 0 ) {
         (void) mach_timebase_info(&sTimebaseInfo);
     }
-    
+
     t2 = mach_absolute_time();
     elapsedNano = (t2-t1) * sTimebaseInfo.numer / sTimebaseInfo.denom;
     printf(str, (float)elapsedNano/1000000.f );

src/owOpenCLSolver.cpp

@@ -40,7 +40,7 @@
 extern int numOfElasticP;
 extern int numOfMembranes;
 
-int myCompare( const void * v1, const void * v2 ); 
+int myCompare( const void * v1, const void * v2 );
 
 /** Constructor of class owOpenCLSolver
  *
@@ -81,7 +81,7 @@ owOpenCLSolver::owOpenCLSolver(const float * position_cpp, const float * velocit
 		create_ocl_kernel("hashParticles", hashParticles);
 		create_ocl_kernel("indexx", indexx);
 		create_ocl_kernel("sortPostPass", sortPostPass);
-		// Additional PCISPH-related kernels 
+		// Additional PCISPH-related kernels
 		create_ocl_kernel("pcisph_computeForcesAndInitPressure", pcisph_computeForcesAndInitPressure);
 		create_ocl_kernel("pcisph_integrate", pcisph_integrate);
 		create_ocl_kernel("pcisph_predictPositions", pcisph_predictPositions);
@@ -215,7 +215,7 @@ void owOpenCLSolver::initializeOpenCL(owConfigProrerty * config)
 		if (ciErrNum == CL_SUCCESS)
 		{
 			printf(" CL_PLATFORM_VERSION [%d]: \t%s\n", i, cBuffer);
-		} 
+		}
 		else
 		{
 			printf(" Error %i in clGetPlatformInfo Call !!!\n\n", ciErrNum);
@@ -286,7 +286,7 @@ void owOpenCLSolver::initializeOpenCL(owConfigProrerty * config)
 	if(result == CL_SUCCESS) std::cout << "CL_CONTEXT_PLATFORM [" << plList <<"]: CL_DEVICE_GLOBAL_MEM_CACHE_SIZE [" << deviceNum <<"]:\t" << val2 <<std::endl;
 	result = devices[deviceNum].getInfo(CL_DEVICE_LOCAL_MEM_SIZE,&val2);
 	if(result == CL_SUCCESS) std::cout << "CL_CONTEXT_PLATFORM " << plList <<": CL_DEVICE_LOCAL_MEM_SIZE ["<< deviceNum <<"]:\t" << val2 << std::endl;
-	
+
 	queue = cl::CommandQueue( context, devices[ deviceNum ], 0, &err );
 	if( err != CL_SUCCESS ){
 		throw std::runtime_error( "failed to create command queue" );
@@ -432,7 +432,7 @@ unsigned int owOpenCLSolver::_runIndexPostPass(owConfigProrerty * config)
 	for(int i=config->gridCellCount;i>=0;i--)
 	{
 		if(gridNextNonEmptyCellBuffer[i] == NO_CELL_ID)
-			gridNextNonEmptyCellBuffer[i] = recentNonEmptyCell; 
+			gridNextNonEmptyCellBuffer[i] = recentNonEmptyCell;
 		else recentNonEmptyCell = gridNextNonEmptyCellBuffer[i];
 	}
 	int err = copy_buffer_to_device( gridNextNonEmptyCellBuffer,gridCellIndexFixedUp,(config->gridCellCount+1) * sizeof( unsigned int ) * 1 );
@@ -532,12 +532,12 @@ unsigned int owOpenCLSolver::_runFindNeighbors(owConfigProrerty * config)
 	int err = queue.enqueueNDRangeKernel(
 		findNeighbors, cl::NullRange, cl::NDRange( (int) (  config->getParticleCount_RoundUp() ) ),
 #if defined( __APPLE__ )
-		cl::NullRange, NULL, NULL );/* 
+		cl::NullRange, NULL, NULL );/*
 		local_work_size can also be a NULL
 		value in which case the OpenCL implementation will
-		determine how to be break the global work-items 
+		determine how to be break the global work-items
 		into appropriate work-group instances.
-		http://www.khronos.org/registry/cl/specs/opencl-1.0.43.pdf, page 109 
+		http://www.khronos.org/registry/cl/specs/opencl-1.0.43.pdf, page 109
 		*/
 #else
 		cl::NDRange( (int)( local_NDRange_size ) ), NULL, NULL );
@@ -918,7 +918,7 @@ unsigned int owOpenCLSolver::_run_computeInteractionWithMembranes_finalize(owCon
  *  @return value taking after enqueue a command to execute a kernel on a device.
  *  More info here (http://www.khronos.org/registry/cl/sdk/1.0/docs/man/xhtml/clEnqueueNDRangeKernel.html)
  */
-unsigned int owOpenCLSolver::_run_pcisph_integrate(int iterationCount, owConfigProrerty * config)
+unsigned int owOpenCLSolver::_run_pcisph_integrate(int iterationCount, int pcisph_integrate_mode, owConfigProrerty * config)
 {
 	// Stage Integrate
 	pcisph_integrate.setArg( 0, acceleration );
@@ -944,6 +944,7 @@ unsigned int owOpenCLSolver::_run_pcisph_integrate(int iterationCount, owConfigP
 	pcisph_integrate.setArg( 20, neighborMap );
 	pcisph_integrate.setArg( 21, config->getParticleCount() );
 	pcisph_integrate.setArg( 22, iterationCount );
+	pcisph_integrate.setArg( 23, pcisph_integrate_mode );
 	int err = queue.enqueueNDRangeKernel(
 		pcisph_integrate, cl::NullRange, cl::NDRange( (int) (  config->getParticleCount_RoundUp() ) ),
 #if defined( __APPLE__ )
@@ -1025,7 +1026,7 @@ void owOpenCLSolver::create_ocl_buffer(const char *name, cl::Buffer &b, const cl
  */
 int owOpenCLSolver::copy_buffer_to_device(const void *host_b, cl::Buffer &ocl_b, const int size )
 {
-	//Actualy we should check  size and type 
+	//Actualy we should check  size and type
 	int err = queue.enqueueWriteBuffer( ocl_b, CL_TRUE, 0, size, host_b );
 	if( err != CL_SUCCESS ){
 		throw std::runtime_error( "Could not enqueue write" );
@@ -1047,7 +1048,7 @@ int owOpenCLSolver::copy_buffer_to_device(const void *host_b, cl::Buffer &ocl_b,
  */
 int owOpenCLSolver::copy_buffer_from_device(void *host_b, const cl::Buffer &ocl_b, const int size )
 {
-	//Actualy we should check  size and type 
+	//Actualy we should check  size and type
 	int err = queue.enqueueReadBuffer( ocl_b, CL_TRUE, 0, size, host_b );
 	if( err != CL_SUCCESS ){
 		throw std::runtime_error( "Could not enqueue read" );

src/owPhysicsFluidSimulator.cpp

@@ -63,6 +63,7 @@ owPhysicsFluidSimulator::owPhysicsFluidSimulator(owHelper * helper,DEVICE dev_ty
 		iterationCount = 0;
 		config = new owConfigProrerty();
 		config->setDeviceType(dev_type);
+		config->integration_method = EULER;//LEAPFROG;
 		// LOAD FROM FILE
 		owHelper::preLoadConfiguration(numOfMembranes, config, numOfLiquidP, numOfElasticP, numOfBoundaryP);
 
@@ -91,8 +92,8 @@ owPhysicsFluidSimulator::owPhysicsFluidSimulator(owHelper * helper,DEVICE dev_ty
 		//The buffers listed below are only for usability and debug
 		density_cpp = new float[ 1 * config->getParticleCount() ];
 		particleIndex_cpp = new unsigned int[config->getParticleCount() * 2];
-		
-		// LOAD FROM FILE	
+
+		// LOAD FROM FILE
 		owHelper::loadConfiguration( position_cpp, velocity_cpp, elasticConnectionsData_cpp, numOfLiquidP, numOfElasticP, numOfBoundaryP, numOfElasticConnections, numOfMembranes,membraneData_cpp, particleMembranesList_cpp, config );		//Load configuration from file to buffer
 
 		if(numOfElasticP != 0){
@@ -185,6 +186,10 @@ 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
+			//Calc next time (t+dt) positions x(t+dt)
+			ocl_solver->_run_pcisph_integrate(iterationCount,0/*=positions_mode*/, config);
+		}
 		ocl_solver->_run_pcisph_computeDensity(config);
 		ocl_solver->_run_pcisph_computeForcesAndInitPressure(config);
 		ocl_solver->_run_pcisph_computeElasticForces(config);
@@ -197,7 +202,13 @@ double owPhysicsFluidSimulator::simulationStep(const bool load_to)
 			iter++;
 		}while( iter < maxIteration );
 
-		ocl_solver->_run_pcisph_integrate(iterationCount, config);		helper->watch_report("_runPCISPH: \t\t%9.3f ms\t3 iteration(s)\n");
+		//and finally calculate v(t+dt)
+		if(config->integration_method == 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{
+			ocl_solver->_run_pcisph_integrate(iterationCount, 2,config);		helper->watch_report("_runPCISPH: \t\t%9.3f ms\t3 iteration(s)\n");
+		}
 		//Handling of Interaction with membranes
 		if(numOfMembranes > 0){
 			ocl_solver->_run_clearMembraneBuffers(config);