diff --git a/inc/owConfigProperty.h b/inc/owConfigProperty.h index 8f1ce64a..8ae24a85 100644 --- a/inc/owConfigProperty.h +++ b/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 }; diff --git a/inc/owOpenCLConstant.h b/inc/owOpenCLConstant.h index 8c7eb563..d45594fb 100644 --- a/inc/owOpenCLConstant.h +++ b/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 diff --git a/inc/owOpenCLSolver.h b/inc/owOpenCLSolver.h index aa877d51..bfdcc3a2 100644 --- a/inc/owOpenCLSolver.h +++ b/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 diff --git a/src/owHelper.cpp b/src/owHelper.cpp index 53030c4a..82e21994 100644 --- a/src/owHelper.cpp +++ b/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 ); diff --git a/src/owOpenCLSolver.cpp b/src/owOpenCLSolver.cpp index ddaf1032..0c440106 100644 --- a/src/owOpenCLSolver.cpp +++ b/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 <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" ); diff --git a/src/owPhysicsFluidSimulator.cpp b/src/owPhysicsFluidSimulator.cpp index b8e28a04..c83b9253 100644 --- a/src/owPhysicsFluidSimulator.cpp +++ b/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); diff --git a/src/owWorldSimulation.cpp b/src/owWorldSimulation.cpp index 3bad43e6..58d53f4f 100644 --- a/src/owWorldSimulation.cpp +++ b/src/owWorldSimulation.cpp @@ -156,7 +156,7 @@ void display(void) helper->refreshTime(); } - glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); + glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT); drawScene(); glPointSize(3.f); glBegin(GL_POINTS); @@ -185,8 +185,8 @@ void display(void) if((int)p_cpp[i*4 + 3] != BOUNDARY_PARTICLE /*&& (int)p_cpp[i*4 + 3] != ELASTIC_PARTICLE*/) { glBegin(GL_POINTS); - if((int)p_cpp[i*4+3]==2) - { + if((int)p_cpp[i*4+3]==2) + { glColor4f( 0, 0, 0, 1.0f);// color of elastic particles glPointSize(6.f); } @@ -220,17 +220,17 @@ void display(void) if((j=(int)ec_cpp[ 4 * i_ec + 0 ])>=0) { i = (i_ec / MAX_NEIGHBOR_COUNT);// + (generateInitialConfiguration!=1)*numOfBoundaryP; - if(i1.f)//muscles + if(ec_cpp[ 4 * i_ec + 2 ]>1.f)//muscles { glLineWidth((GLfloat)1.0); - if(ec_cpp[4*i_ec+2]-floor(ec_cpp[4*i_ec+2])>0.45f) + if(ec_cpp[4*i_ec+2]-floor(ec_cpp[4*i_ec+2])>0.45f) { if(muscle_activation_signal_cpp[ (int)(floor( ec_cpp[4*i_ec+2])-1) ]>0.1) glLineWidth((GLfloat)6.0); else glLineWidth((GLfloat)2.0); - glColor4b(127/2, 0, 255/2, 255/2);/* muscle_number+0.5 <--> violet*/ + glColor4b(127/2, 0, 255/2, 255/2);/* muscle_number+0.5 <--> violet*/ glBegin(GL_LINES); glVertex3f( (p_cpp[i*4+0]-loacalConfig->xmax/2)*sc , (p_cpp[i*4+1]-loacalConfig->ymax/2)*sc, (p_cpp[i*4+2]-loacalConfig->zmax/2)*sc ); glColor4b(255/2, 255/2, 255/2, 255/2); @@ -238,11 +238,11 @@ void display(void) glEnd(); } else - if(ec_cpp[4*i_ec+2]-floor(ec_cpp[4*i_ec+2])>0.35f) - { + if(ec_cpp[4*i_ec+2]-floor(ec_cpp[4*i_ec+2])>0.35f) + { if(muscle_activation_signal_cpp[ (int)(floor( ec_cpp[4*i_ec+2])-1) ]>0.1) glLineWidth((GLfloat)6.0); else glLineWidth((GLfloat)2.0); - glColor4b(255/2, 0, 255/2, 255/2);/* muscle_number+0.4 <--> magenta*/ + glColor4b(255/2, 0, 255/2, 255/2);/* muscle_number+0.4 <--> magenta*/ glBegin(GL_LINES); glVertex3f( (p_cpp[i*4+0]-loacalConfig->xmax/2)*sc , (p_cpp[i*4+1]-loacalConfig->ymax/2)*sc, (p_cpp[i*4+2]-loacalConfig->zmax/2)*sc ); glColor4b(255/2, 255/2, 255/2, 255/2); @@ -250,11 +250,11 @@ void display(void) glEnd(); } else - if(ec_cpp[4*i_ec+2]-floor(ec_cpp[4*i_ec+2])>0.25f) - { + if(ec_cpp[4*i_ec+2]-floor(ec_cpp[4*i_ec+2])>0.25f) + { if(muscle_activation_signal_cpp[ (int)(floor( ec_cpp[4*i_ec+2])-1) ]>0.1) glLineWidth((GLfloat)6.0); else glLineWidth((GLfloat)2.0); - glColor4b(255/2, 127/2, 0, 255/2);/* muscle_number+0.3 <--> orange*/ + glColor4b(255/2, 127/2, 0, 255/2);/* muscle_number+0.3 <--> orange*/ glBegin(GL_LINES); glVertex3f( (p_cpp[i*4+0]-loacalConfig->xmax/2)*sc , (p_cpp[i*4+1]-loacalConfig->ymax/2)*sc, (p_cpp[i*4+2]-loacalConfig->zmax/2)*sc ); glColor4b(255/2, 255/2, 255/2, 255/2); @@ -262,11 +262,11 @@ void display(void) glEnd(); } else - if(ec_cpp[4*i_ec+2]-floor(ec_cpp[4*i_ec+2])>0.15f) - { + if(ec_cpp[4*i_ec+2]-floor(ec_cpp[4*i_ec+2])>0.15f) + { if(muscle_activation_signal_cpp[ (int)(floor( ec_cpp[4*i_ec+2])-1) ]>0.1) glLineWidth((GLfloat)6.0); else glLineWidth((GLfloat)2.0); - glColor4b(255/2, 0, 0, 255/2);/* muscle_number+0.2 <--> red*/ + glColor4b(255/2, 0, 0, 255/2);/* muscle_number+0.2 <--> red*/ glBegin(GL_LINES); glVertex3f( (p_cpp[i*4+0]-loacalConfig->xmax/2)*sc , (p_cpp[i*4+1]-loacalConfig->ymax/2)*sc, (p_cpp[i*4+2]-loacalConfig->zmax/2)*sc ); glColor4b(255/2, 255/2, 255/2, 255/2); @@ -295,7 +295,7 @@ void display(void) glVertex3f( (p_cpp[j*4+0]-loacalConfig->xmax/2)*sc , (p_cpp[j*4+1]-loacalConfig->ymax/2)*sc, (p_cpp[j*4+2]-loacalConfig->zmax/2)*sc ); glEnd(); } - + ecc++; } } @@ -333,11 +333,11 @@ void display(void) inline void drawScene() { // [7]----[6] - // / | /| - // [3]----[2] | - // | [4]--|-[5] + // / | /| + // [3]----[2] | + // | [4]--|-[5] // | / | / - // [0]----[1] + // [0]----[1] Vector3D vcenter(0,0,0); Vector3D vbox[8]; float s_v = 1 /(simulationScale);// = 1 m in simulation @@ -361,7 +361,7 @@ inline void drawScene() vbox[7] = Vector3D(loacalConfig->xmin,loacalConfig->ymax,loacalConfig->zmax); // Display user interface if enabled bool displayInfos = true; - if (displayInfos) + if (displayInfos) { glDisable(GL_DEPTH_TEST); glBlendFunc(GL_ONE_MINUS_DST_COLOR, GL_ZERO); // invert color @@ -394,9 +394,9 @@ inline void drawScene() v7 = Vector3D( loacalConfig->xmax/2, loacalConfig->ymax/2, loacalConfig->zmax/2)*sc; v8 = Vector3D( -loacalConfig->xmax/2, loacalConfig->ymax/2, loacalConfig->zmax/2)*sc; glColor3ub(255,255,255);//yellow - glVertex3d(v1.x,v1.y,v1.z); + glVertex3d(v1.x,v1.y,v1.z); glVertex3d(v2.x,v2.y,v2.z); - + glColor3ub(255,255,255);//yellow glVertex3d(v2.x,v2.y,v2.z); @@ -436,7 +436,7 @@ inline void drawScene() // glBegin(GL_LINES); glColor3ub(0,0,0);//black - + Vector3D v_s = Vector3D( -loacalConfig->xmax/2 + s_v, loacalConfig->ymax/2, loacalConfig->zmax/2)*sc; glVertex3d(v_s.x, v_s.y, v_s.z); @@ -482,7 +482,7 @@ void renderInfo(int x, int y) numOfElasticP, numOfBoundaryP,loacalConfig->getParticleCount()); glPrint( 0 , 2 , label, m_font); - glColor3f (1.0F, 1.0F, 1.0F); + glColor3f (1.0F, 1.0F, 1.0F); sprintf(label,"Selected device: %s FPS = %.2f, time step: %d (%f s)", device_full_name+7, fps, fluid_simulation->getIteration(),((float)fluid_simulation->getIteration())*timeStep); glPrint( 0 , 17 , label, m_font); } @@ -501,7 +501,7 @@ void renderInfo(int x, int y) glVertex2f((GLfloat) s_v,(GLfloat)y_m + 5.f ); glEnd(); glPrint( s_v , y_m + 15.f , "1E-02 m", m_font); - glBegin(GL_LINES); + glBegin(GL_LINES); glVertex2f((GLfloat) s_v_10,(GLfloat)y_m + 0.f); glVertex2f((GLfloat) s_v_10,(GLfloat)y_m + 5.f); glEnd(); @@ -516,7 +516,7 @@ void renderInfo(int x, int y) } if(flag){ for(int i = 1;i <= count_s; i++){ - glBegin(GL_LINES); + glBegin(GL_LINES); glVertex2f((GLfloat) s_v/pow(10.f,i + 1),(GLfloat)y_m + 0.f); glVertex2f((GLfloat) s_v/pow(10.f,i + 1),(GLfloat)y_m + 5.f); glEnd(); @@ -552,10 +552,10 @@ void respond_mouse(int button, int state, int x, int y) buttonState = 3; int mods; mods = glutGetModifiers(); - if (mods & GLUT_ACTIVE_CTRL) + if (mods & GLUT_ACTIVE_CTRL) { buttonState = 2; - } + } if(state == GLUT_UP) buttonState = 0; old_x=x; @@ -576,12 +576,12 @@ void respond_mouse(int button, int state, int x, int y) // GLUT callback // called on mouse movement -void mouse_motion (int x, int y) +void mouse_motion (int x, int y) { float dx,dy; - dy = (float)(y - old_y); + dy = (float)(y - old_y); dx = (float)(x - old_x); - + if(buttonState == 1) { camera_rot[0] += dy / 5.0f; @@ -663,11 +663,11 @@ void RespondKey(unsigned char key, int x, int y) * animation, this idle() function must update not only the * main window but also all derived subwindows */ -void idle (void) -{ - glutSetWindow (winIdMain); - glutPostRedisplay (); -} +void idle (void) +{ + glutSetWindow (winIdMain); + glutPostRedisplay (); +} //static char label[1000]; /* Storage for current string */ void Timer(int value) @@ -679,13 +679,13 @@ void Timer(int value) GLvoid resize(GLsizei width, GLsizei height){ - if(height == 0) { height = 1; } - if(width == 0) { width = 1; } - + if(height == 0) { height = 1; } + if(width == 0) { width = 1; } + glViewport(0, 0, width, height); // Set view area glMatrixMode(GL_PROJECTION); glLoadIdentity(); - + float aspectRatio = (GLfloat)width / (GLfloat)height; if (aspectRatio>1.f) glFrustum(-1*aspectRatio, 1*aspectRatio, -1, 1, 3, 45); @@ -743,14 +743,14 @@ void run(int argc, char** argv, const bool with_graphics, const bool load_to) helper = new owHelper(); if(!load_from_file){ DEVICE dev_type = CPU; - for(int i = 1; igetConfig(); + loacalConfig = fluid_simulation->getConfig(); } else{ loacalConfig = new owConfigProrerty(); @@ -767,7 +767,7 @@ void run(int argc, char** argv, const bool with_graphics, const bool load_to) glutInitWindowSize(1024, 1024); glutInitWindowPosition(100, 100); winIdMain = glutCreateWindow("Palyanov Andrey for OpenWorm: OpenCL PCISPH fluid + elastic matter + membranes [2013]: C.elegans body generator demo"); - glutIdleFunc (idle); + glutIdleFunc (idle); //Init physic Simulation init(); glutDisplayFunc(display); diff --git a/src/sphFluid.cl b/src/sphFluid.cl index d3c89a06..4cef99c0 100644 --- a/src/sphFluid.cl +++ b/src/sphFluid.cl @@ -40,8 +40,8 @@ * Boundary handling - [4] M. Ihmsen, N. Akinci, M. Gissler, M. Teschner, * Boundary Handling and Adaptive Time-stepping for PCISPH * Proc. VRIPHYS, Copenhagen, Denmark, pp. 79-88, Nov 11-12, 2010 - * Surface Tension - [5] M. Becker, M. Teschner. Weakly compressible SPH for free surface flows - * // Proceedings of the 2007 ACM SIGGRAPH/Eurographics + * Surface Tension - [5] M. Becker, M. Teschner. Weakly compressible SPH for free surface flows + * // Proceedings of the 2007 ACM SIGGRAPH/Eurographics * symposium on Computer animation, pages 209-217. */ @@ -79,7 +79,7 @@ #endif #ifdef cl_amd_printf #pragma OPENCL EXTENSION cl_amd_printf : enable -//#define PRINTF_ON // this comment because using printf leads to very slow work on Radeon r9 290x on my machine +//#define PRINTF_ON // this comment because using printf leads to very slow work on Radeon r9 290x on my machine // don't know why #elif defined(cl_intel_printf) #pragma OPENCL EXTENSION cl_intel_printf : enable @@ -125,7 +125,7 @@ __kernel void clearBuffers( nm[ outIdx++ ] = fdata; } -int cellId( +int cellId( int4 cellFactors_, int gridCellsX, int gridCellsY, @@ -147,7 +147,7 @@ int cellId( float hashGridCellSizeInv ) { - //xmin, ymin, zmin + //xmin, ymin, zmin int4 result; result.x = (int)( position.x * hashGridCellSizeInv ); result.y = (int)( position.y * hashGridCellSizeInv ); @@ -170,7 +170,7 @@ __kernel void hashParticles( int id = get_global_id( 0 ); if( id >= PARTICLE_COUNT ) return; float4 _position = position[ id ]; - int4 cellFactors_ = cellFactors( _position, xmin, ymin, zmin, hashGridCellSizeInv ); + int4 cellFactors_ = cellFactors( _position, xmin, ymin, zmin, hashGridCellSizeInv ); int cellId_ = cellId( cellFactors_, gridCellsX, gridCellsY, gridCellsZ ) & 0xffff; // truncate to low 16 bits uint2 result; PI_CELL_ID( result ) = cellId_; @@ -200,7 +200,7 @@ __kernel void sortPostPass( { int id = get_global_id( 0 ); if( id >= PARTICLE_COUNT ) return; - uint2 spi = particleIndex[ id ];//contains id of cell and id of particle it has sorted + uint2 spi = particleIndex[ id ];//contains id of cell and id of particle it has sorted int serialId = PI_SERIAL_ID( spi );//get a particle Index int cellId = PI_CELL_ID( spi );//get a cell Index float4 position_ = position[ serialId ];//get position by serialId @@ -211,7 +211,7 @@ __kernel void sortPostPass( particleIndexBack[ serialId ] = id; } /** Calculating start position in particleIndex for every cell - * Kernel fill up gridCellIndex buffer empty cell + * Kernel fill up gridCellIndex buffer empty cell * (spatial cell which has no particle inside at this time is filling by -1). */ __kernel void indexx( @@ -230,7 +230,7 @@ __kernel void indexx( // add the nth+1 index value gridCellIndex[ id ] = PARTICLE_COUNT; return; - } + } if( id == 0 ){ gridCellIndex[ id ] = 0; return; @@ -255,7 +255,7 @@ __kernel void indexx( bool isLow = ( sampleCellId < id ); low = isLow ? idx + 1 : low; bool isMiddle = !( isHigh || isLow ); - + bool zeroCase = ( idx == 0 && isMiddle ); int sampleM1CellId = zeroCase ? -1 : PI_CELL_ID( sampleMinus1 ); converged = isMiddle && ( zeroCase || sampleM1CellId < sampleCellId ); @@ -280,7 +280,7 @@ int getMaxIndex( } // Neighbour Search algorithm function and kernel block /** Searching for neighbour in particular spatial cell for particular particle - * It takes every particles from particular cell and check if it satisfy + * It takes every particles from particular cell and check if it satisfy * a condition that distance between particles is <= closest_distance */ int searchForNeighbors_b( @@ -296,7 +296,7 @@ int searchForNeighbors_b( int *found_count ) { - int baseParticleId = gridCellIndex[ searchCell_ ]; + int baseParticleId = gridCellIndex[ searchCell_ ]; int nextParticleId = gridCellIndex[ searchCell_ + 1 ]; int particleCountThisCell = nextParticleId - baseParticleId; // Calcuating how many particle is containining in particular cell int i = 0; @@ -328,13 +328,13 @@ int searchForNeighbors_b( /** Return value of cellId from gridCellIndexFixedUp * for particular cell and offset (deltaX,Y,Z) */ -int searchCell( +int searchCell( int cellId, int deltaX, int deltaY, int deltaZ, - int gridCellsX, - int gridCellsY, + int gridCellsX, + int gridCellsY, int gridCellsZ, int gridCellCount ) @@ -384,72 +384,72 @@ __kernel void findNeighbors( closest_distances[k] = r_thr2; closest_indexes[k] = -1; } - + searchCells[0] = myCellId; - + // p is the current particle position within the bounds of the hash grid float4 p; float4 p0 = (float4)( xmin, ymin, zmin, 0.0f ); p = position_ - p0; - + // cf is the min,min,min corner of the current cell int4 cellFactors_ = cellFactors( position_, xmin, ymin, zmin, hashGridCellSizeInv ); float4 cf; cf.x = cellFactors_.x * hashGridCellSize; cf.y = cellFactors_.y * hashGridCellSize; cf.z = cellFactors_.z * hashGridCellSize; - + // lo.A is true if the current position is in the low half of the cell for dimension A int4 lo; lo = (( p - cf ) < h ); - + int4 delta; int4 one = (int4)( 1, 1, 1, 1 ); delta = one + 2 * lo; - // search surrounding cells 1..8 + // search surrounding cells 1..8 searchCells[1] = searchCell( myCellId, delta.x, 0, 0, gridCellsX, gridCellsY, gridCellsZ, gridCellCount ); searchCells[2] = searchCell( myCellId, 0, delta.y, 0, gridCellsX, gridCellsY, gridCellsZ, gridCellCount ); searchCells[3] = searchCell( myCellId, 0, 0, delta.z, gridCellsX, gridCellsY, gridCellsZ, gridCellCount ); searchCells[4] = searchCell( myCellId, delta.x, delta.y, 0, gridCellsX, gridCellsY, gridCellsZ, gridCellCount ); searchCells[5] = searchCell( myCellId, delta.x, 0, delta.z, gridCellsX, gridCellsY, gridCellsZ, gridCellCount ); searchCells[6] = searchCell( myCellId, 0, delta.y, delta.z, gridCellsX, gridCellsY, gridCellsZ, gridCellCount ); - searchCells[7] = searchCell( myCellId, delta.x, delta.y, delta.z, gridCellsX, gridCellsY, gridCellsZ, gridCellCount ); - + searchCells[7] = searchCell( myCellId, delta.x, delta.y, delta.z, gridCellsX, gridCellsY, gridCellsZ, gridCellCount ); + int last_farthest = 0; // Search neighbour particles in every cells from searchCells list last_farthest = searchForNeighbors_b( searchCells[0], gridCellIndex, position_, id, sortedPosition, neighborMap, closest_indexes, closest_distances, last_farthest, &found_count ); - + last_farthest = searchForNeighbors_b( searchCells[1], gridCellIndex, position_, id, sortedPosition, neighborMap, closest_indexes, closest_distances, last_farthest, &found_count ); - + last_farthest = searchForNeighbors_b( searchCells[2], gridCellIndex, position_, id, sortedPosition, neighborMap, closest_indexes, closest_distances, last_farthest, &found_count ); - + last_farthest = searchForNeighbors_b( searchCells[3], gridCellIndex, position_, id, sortedPosition, neighborMap, closest_indexes, closest_distances, last_farthest, &found_count ); - + last_farthest = searchForNeighbors_b( searchCells[4], gridCellIndex, position_, id, sortedPosition, neighborMap, closest_indexes, closest_distances, last_farthest, &found_count ); - + last_farthest = searchForNeighbors_b( searchCells[5], gridCellIndex, position_, id, sortedPosition, neighborMap, closest_indexes, closest_distances, last_farthest, &found_count ); - + last_farthest = searchForNeighbors_b( searchCells[6], gridCellIndex, position_, id, sortedPosition, neighborMap, closest_indexes, closest_distances, last_farthest, &found_count ); - + last_farthest = searchForNeighbors_b( searchCells[7], gridCellIndex, position_, id, sortedPosition, neighborMap, closest_indexes, closest_distances, last_farthest, &found_count ); - // storing all found neighbors into neighborMap buffer + // storing all found neighbors into neighborMap buffer for(int j=0; j=PARTICLE_COUNT) return; position[PARTICLE_COUNT + id] = (float4)(0.0f,0.0f,0.0f,0.0f); //extra memory to store changes in considered particles's coordinates due to interaction with membranes. Need to make it zero every step. @@ -986,7 +986,7 @@ __kernel void clearMembraneBuffers( */ float calcDeterminant3x3(float4 c1, float4 c2, float4 c3) {// here we expect the following structure of input vectors (0-th component of each is not used, 1,2,3 - used) - return c1.x*c2.y*c3.z + c1.y*c2.z*c3.x + c1.z*c2.x*c3.y + return c1.x*c2.y*c3.z + c1.y*c2.z*c3.x + c1.z*c2.x*c3.y - c1.z*c2.y*c3.x - c1.x*c2.z*c3.y - c1.y*c2.x*c3.z; } @@ -1030,13 +1030,13 @@ float4 calculateProjectionOfPointToPlane(float4 ps, float4 pa, float4 pb, float4 #ifdef PRINTF_ON printf("\ndenominator equal to zero\n"); #endif - pm.w = -1;//indicates error + pm.w = -1;//indicates error } return pm; } float calculateTriangleSquare(float4 v1, float4 v2, float4 v3) { - // here 'v' is for vertex or vector, anyway v1, v2, v3 are coordinates of 3 points in 3D. + // here 'v' is for vertex or vector, anyway v1, v2, v3 are coordinates of 3 points in 3D. // 4-th coordinate is not used // first calc 2 vectors: v21 and v31 float4 a = v2 - v1;//v21 @@ -1060,16 +1060,16 @@ __kernel void computeInteractionWithMembranes( int numOfElasticP, float r0 ) { - int id = get_global_id( 0 ); + int id = get_global_id( 0 ); if(id>=PARTICLE_COUNT) return; - id = particleIndexBack[id]; + id = particleIndexBack[id]; int id_source_particle = PI_SERIAL_ID( particleIndex[id] ); int jd_source_particle; //float4 position_ = sortedPosition[ id ]; float4 position_ = position[ id ]; if((int)(position[ id_source_particle ].w) == BOUNDARY_PARTICLE) return; - if((int)(position[ id_source_particle ].w) != LIQUID_PARTICLE) return; //!!! currently we consider only liquid particles - //!!! interacting with membranes + if((int)(position[ id_source_particle ].w) != LIQUID_PARTICLE) return; //!!! currently we consider only liquid particles + //!!! interacting with membranes int jd, idx = id * MAX_NEIGHBOR_COUNT; int mdi;//membraneData index int i,j,k;//these i and j have nothing in common with id and jd indexes @@ -1106,7 +1106,7 @@ __kernel void computeInteractionWithMembranes( vector_id_jd = position[id_source_particle] - position[jd_source_particle]; vector_id_jd.z = 0.0f; //mv change from subscripting _distance_id_jd = sqrt(dot(vector_id_jd,vector_id_jd)); - // elastic matter particles have no information + // elastic matter particles have no information // about participation in membrane composition // Let's get it - check corresponding position of particleMembranesList (if it is non-empty) for(int mli=0;mli0.0f) { normal_to_ijk_plane /= normal_to_ijk_plane_length;// normalized now membrane_jd_normal_vector[membrane_jd_counter] += normal_to_ijk_plane; membrane_ijk_counter++; - // so, we consider i-th particle and a number of its neighbours which belong to membrane(s). + // so, we consider i-th particle and a number of its neighbours which belong to membrane(s). // normal vectors are calculated for all of them. // now it's time to calculate forces: // 1) force F_i, acting on i-th particle @@ -1157,8 +1157,8 @@ __kernel void computeInteractionWithMembranes( } // ok, we finally have projection of considered particle on the plane of i-j-k triangle. // If triangle's square >0 and if projection point is inside the triangle (not outside) - // then this triangle is located is such way that we have to take it into account and - // calculate repulsion from it. + // then this triangle is located is such way that we have to take it into account and + // calculate repulsion from it. } else break; }//22222222222222222222222222222222222222222222222222222222222222222222222222222222222 @@ -1177,11 +1177,11 @@ __kernel void computeInteractionWithMembranes( if(membrane_jd_counter>0) { int nc = 0; - float4 n_c_i = (float4)(0.f,0.f,0.f,0.f); + float4 n_c_i = (float4)(0.f,0.f,0.f,0.f); float4 n_m; float w_c_im, w_c_im_sum = 0.f, w_c_im_second_sum = 0.f; float4 delta_pos; - float4 velocity_membrane_average = (float4)(0.f,0.f,0.f,0.f); + float4 velocity_membrane_average = (float4)(0.f,0.f,0.f,0.f); float n_c_i_length,x_im_dist; int id_m_source_particle;//index of i-th particle's (current) neighbours which are membrane particles // they are already in the list @@ -1230,18 +1230,18 @@ __kernel void computeInteractionWithMembranes_finalize( int PARTICLE_COUNT ) { - int id = get_global_id( 0 ); + int id = get_global_id( 0 ); if(id>=PARTICLE_COUNT) return; - - id = particleIndexBack[id]; + + id = particleIndexBack[id]; int id_source_particle = PI_SERIAL_ID( particleIndex[id] ); int jd_source_particle; float4 position_ = position[ id ]; float v2; if((int)(position[ id_source_particle ].w) == BOUNDARY_PARTICLE) return; - //if((int)(position[ id_source_particle ].w) != LIQUID_PARTICLE) return; //!!! currently we consider only liquid particles - //!!! interacting with membranes + //if((int)(position[ id_source_particle ].w) != LIQUID_PARTICLE) return; //!!! currently we consider only liquid particles + //!!! interacting with membranes position[ id_source_particle ] += position[ PARTICLE_COUNT + id_source_particle ]; /* velocity[ PARTICLE_COUNT + id_source_particle ].w = 0; @@ -1279,29 +1279,38 @@ __kernel void pcisph_integrate( float r0, __global float2 * neighborMap, int PARTICLE_COUNT, - int iterationCount + int iterationCount, + int mode ) { - int id = get_global_id( 0 ); + int id = get_global_id( 0 ); if(id>=PARTICLE_COUNT) return; - id = particleIndexBack[id]; + id = particleIndexBack[id]; int id_source_particle = PI_SERIAL_ID( particleIndex[id] ); if((int)(position[ id_source_particle ].w) == BOUNDARY_PARTICLE) { return; } - float4 acceleration_t = acceleration[ PARTICLE_COUNT*2+id_source_particle ]; acceleration_t.w = 0.f; - float4 acceleration_t_dt = acceleration[ id ] + acceleration[ PARTICLE_COUNT+id ]; acceleration_t_dt.w = 0.f; + //float4 acceleration_t = acceleration[ PARTICLE_COUNT*2+id_source_particle ]; acceleration_t.w = 0.f; + //float4 acceleration_t_dt = acceleration[ id ] + acceleration[ PARTICLE_COUNT+id ]; acceleration_t_dt.w = 0.f; + //float4 velocity_t = sortedVelocity[ id ]; + //float4 position_t = sortedPosition[ id ]; + if(iterationCount==0) + { + acceleration[ PARTICLE_COUNT*2+id_source_particle ] = + acceleration[ id ] + acceleration[ PARTICLE_COUNT+id ]; + return; + } + float4 acceleration_t = acceleration[ PARTICLE_COUNT*2+id_source_particle ]; acceleration_t.w = 0.f; float4 velocity_t = sortedVelocity[ id ]; - float4 position_t = sortedPosition[ id ]; - if(iterationCount==0) - acceleration_t = (float4)(0.0f,-9.8f,0.0f,0.0f); + float particleType = position[ id_source_particle ].w; + ////////////////////////////////////////////////////////////// -//!!/// so-called "Velocity Verlet" integration - /// (similar to leapfrog method, except that the velocity and position are calculated - /// at the same value of the time variable (Leapfrog does not, as the name suggests). +//!!/// so-called "Velocity Verlet" integration + /// (similar to leapfrog method, except that the velocity and position are calculated + /// at the same value of the time variable (Leapfrog does not, as the name suggests). ///========================================================== - /// 2-nd order of precision for the position: + /// 2-nd order of precision for the position: /// [ O(delta_t^2) global (cumulative) error over a constant interval of time ] /// self-starting, minimizes roundoff errors /// http://www.saylor.org/site/wp-content/uploads/2011/06/MA221-6.1.pdf, page 3 @@ -1319,21 +1328,40 @@ __kernel void pcisph_integrate( //http://www.richardlord.net/presentations/physics-for-flash-games //float4 position_t_dt = position_t + (velocity_t + (velocity_t+acceleration_t*timeStep) )*0.5f*timeStep*simulationScaleInv; //float4 velocity_t_dt = velocity_t + (acceleration_t+acceleration_t_dt)*0.5f*timeStep; - + ////////////////////////////////////////////////////////////////////////////////////////////////// // The semi-implicit Euler is a first-order integrator, just as the standard Euler method. // // This means that it commits a global error of the order of dt. However, the semi-implicit // // Euler method is a symplectic integrator, unlike the standard method. As a consequence, // // the semi-implicit Euler method almost conserves the energy (when the Hamiltonian is // - // time-independent). // + // time-independent). // // http://en.wikipedia.org/wiki/Semi-implicit_Euler // ////////////////////////////////////////////////////////////////////////////////////////////////// // Semi-)Implicit Euler integrator. 1st order, symplectic (has a sort of "global" stability) // // Most of the usual numerical methods, like the primitive Euler scheme // // and the classical Runge-Kutta scheme, are not symplectic integrators. // //!!////////////////////////////////////////////////////////////////////////////////////////////////// - /**/ float4 velocity_t_dt = velocity_t + (acceleration_t_dt)*timeStep; // - /**/ float4 position_t_dt = position_t + (velocity_t_dt)*timeStep*simulationScaleInv; // + if(mode == 2){ + float4 acceleration_t_dt = acceleration[ id ] + acceleration[ PARTICLE_COUNT+id ]; acceleration_t_dt.w = 0.f; + float4 position_t = sortedPosition[ id ]; + float4 velocity_t_dt = velocity_t + (acceleration_t_dt)*timeStep; // + float4 position_t_dt = position_t + (velocity_t_dt)*timeStep*simulationScaleInv; // + + float particleType = position[ id_source_particle ].w; + computeInteractionWithBoundaryParticles(id,r0,neighborMap,particleIndexBack,particleIndex,position,velocity,&position_t_dt, true, &velocity_t_dt,PARTICLE_COUNT); + + + velocity[ id_source_particle ] = velocity_t_dt; + position[ id_source_particle ] = position_t_dt; + position[ id_source_particle ].w = particleType; + //velocity[ id_source_particle ] = (float4)((float)velocity_t_dt_x, (float)velocity_t_dt_y, (float)velocity_t_dt_z, 0.f); + //position[ id_source_particle ] = (float4)((float)position_t_dt_x, (float)position_t_dt_y, (float)position_t_dt_z, particleType); + + acceleration[PARTICLE_COUNT*2+id_source_particle] = acceleration_t_dt; + return; + } + /**/// float4 velocity_t_dt = velocity_t + (acceleration_t_dt)*timeStep; // + /**/// float4 position_t_dt = position_t + (velocity_t_dt)*timeStep*simulationScaleInv; // ////////////////////////////////////////////////////////////////////////////////////////////////// //printf("\n===[ timeStep= %5e ]===",timeStep); ////////////////////////////////////////////////////////////////////////////////////////////////// @@ -1345,12 +1373,44 @@ __kernel void pcisph_integrate( //Find the new momentum based on the force and HALF of the small time step interval (not the whole time step) //Find the new position. //Find the next new momentum with the other half of the time step. - //A second way to write the leapfrog looks quite different at first sight. Defining all quantities only at integer times, we can write: - /**/// float4 position_t_dt = position_t + (velocity_t*timeStep + acceleration_t*timeStep*timeStep/2.f)*simulationScaleInv; - /**/// float4 velocity_t_dt = velocity_t + (acceleration_t + acceleration_t_dt)*timeStep/2.f; + //A second way to write the leapfrog looks quite different at first sight. Defining all quantities only at integer times, we can write: + /**/// float4 position_t_dt = position_t + (velocity_t*timeStep + acceleration_t*timeStep*timeStep/2.f)*simulationScaleInv; + /**/// float4 velocity_t_dt = velocity_t + (acceleration_t + acceleration_t_dt)*timeStep/2.f; // for floats it works with a significant error, which neglects all advantages of this really nice method // for example, at first time step we get 2.17819E-03 instead of 2.18000E-03, and such things occur at every step and accumulate. // switching to doubles. + if(mode==0/*positions_mode*/) + { + float4 position_t = sortedPosition[ id ]; + float4 position_t_dt = position_t + (velocity_t*timeStep + acceleration_t*timeStep*timeStep/2.f)*simulationScaleInv; + sortedPosition[ id ] = position_t_dt; + sortedPosition[ id ].w = particleType; + //printf(" ==> mode0\n"); + //printf(" ==> x(t) %f\n",position_t.y); + //printf(" ==> v(t) %f\n",velocity_t.y); + //printf(" ==> a(t) %f\n",acceleration_t.y); + } + else + if(mode==1/*velocities_mode*/) + { + float4 position_t_dt = sortedPosition[ id ];//necessary for computeInteractionsWithBoundaryParticles() + float4 acceleration_t_dt = acceleration[ id ] + acceleration[ PARTICLE_COUNT+id ]; acceleration_t_dt.w = 0.f; + float4 velocity_t_dt = velocity_t + (acceleration_t + acceleration_t_dt)*timeStep/2.f; + + computeInteractionWithBoundaryParticles(id,r0,neighborMap,particleIndexBack,particleIndex,position,velocity,&position_t_dt, true, &velocity_t_dt,PARTICLE_COUNT); + velocity[ id_source_particle ] = velocity_t_dt; + acceleration[PARTICLE_COUNT*2+id_source_particle] = acceleration_t_dt; + + position[ id_source_particle ] = position_t_dt; + position[ id_source_particle ].w = particleType; + + //printf(" ==> mode1\n"); + //printf(" ==> x(t+dt) %f\n",position_t_dt.y); + //printf(" ==> v(t+dt) %f\n",velocity_t_dt.y); + //printf(" ==> a(t+dt) %f\n",acceleration_t_dt.y); + + } + return; //!!//Explicit Euler (1st order) // float4 velocity_t_dt = velocity_t + (acceleration_t)*timeStep; // float4 position_t_dt = position_t + (velocity_t_dt+velocity_t)*0.5*timeStep*simulationScaleInv; @@ -1369,17 +1429,17 @@ __kernel void pcisph_integrate( float4 position_t_dt = position_t + ( k1_p + k2_p ) * 1.0f/2.0f; */ // in Chao Fang realization here is also acceleration 'speed limit' applied - if(position_t_dt.xxmax-0.000001f) position_t_dt.x = xmax-0.000001f;//A.Palyanov 30.08.2012 if(position_t_dt.y>ymax-0.000001f) position_t_dt.y = ymax-0.000001f;//A.Palyanov 30.08.2012 - if(position_t_dt.z>zmax-0.000001f) position_t_dt.z = zmax-0.000001f;//A.Palyanov 30.08.2012 - // better replace 0.0000001 with smoothingRadius*0.001 or smth like this - float particleType = position[ id_source_particle ].w; + if(position_t_dt.z>zmax-0.000001f) position_t_dt.z = zmax-0.000001f;//A.Palyanov 30.08.2012*/ + // better replace 0.0000001 with smoothingRadius*0.001 or smth like this + /*float particleType = position[ id_source_particle ].w; computeInteractionWithBoundaryParticles(id,r0,neighborMap,particleIndexBack,particleIndex,position,velocity,&position_t_dt, true, &velocity_t_dt,PARTICLE_COUNT); velocity[ id_source_particle ] = velocity_t_dt; position[ id_source_particle ] = position_t_dt; position[ id_source_particle ].w = particleType; - acceleration[PARTICLE_COUNT*2+id_source_particle] = acceleration_t_dt; + acceleration[PARTICLE_COUNT*2+id_source_particle] = acceleration_t_dt;*/ }