Extras/sph/fluids/fluid_system.cpp

   1 /*
   2   FLUIDS v.1 - SPH Fluid Simulator for CPU and GPU
   3   Copyright (C) 2008. Rama Hoetzlein, http://www.rchoetzlein.com
   4
   5   ZLib license
   6   This software is provided 'as-is', without any express or implied
   7   warranty.  In no event will the authors be held liable for any damages
   8   arising from the use of this software.
   9
  10   Permission is granted to anyone to use this software for any purpose,
  11   including commercial applications, and to alter it and redistribute it
  12   freely, subject to the following restrictions:
  13
  14   1. The origin of this software must not be misrepresented; you must not
  15      claim that you wrote the original software. If you use this software
  16      in a product, an acknowledgment in the product documentation would be
  17      appreciated but is not required.
  18   2. Altered source versions must be plainly marked as such, and must not be
  19      misrepresented as being the original software.
  20   3. This notice may not be removed or altered from any source distribution.
  21 */
  22
  23
  24
  25 #include <conio.h>
  26
  27 #ifdef _MSC_VER
  28         #include <gl/glut.h>
  29 #else
  30         #include <GL/glut.h>
  31 #endif
  32
  33 #include "common_defs.h"
  34 #include "mtime.h"
  35 #include "fluid_system.h"
  36
  37 #ifdef BUILD_CUDA
  38         #include "fluid_system_host.cuh"
  39 #endif
  40
  41 #define EPSILON                 0.00001f                        //for collision detection
  42
  43 FluidSystem::FluidSystem ()
  44 {
  45 }
  46
  47 void FluidSystem::Initialize ( int mode, int total )
  48 {
  49         if ( mode != BFLUID ) {
  50                 printf ( "ERROR: FluidSystem not initialized as BFLUID.\n");
  51         }
  52         PointSet::Initialize ( mode, total );
  53
  54         FreeBuffers ();
  55         AddBuffer ( BFLUID, sizeof ( Fluid ), total );
  56         AddAttribute ( 0, "pos", sizeof ( Vector3DF ), false );
  57         AddAttribute ( 0, "color", sizeof ( DWORD ), false );
  58         AddAttribute ( 0, "vel", sizeof ( Vector3DF ), false );
  59         AddAttribute ( 0, "ndx", sizeof ( unsigned short ), false );
  60         AddAttribute ( 0, "age", sizeof ( unsigned short ), false );
  61
  62         AddAttribute ( 0, "pressure", sizeof ( double ), false );
  63         AddAttribute ( 0, "density", sizeof ( double ), false );
  64         AddAttribute ( 0, "sph_force", sizeof ( Vector3DF ), false );
  65         AddAttribute ( 0, "next", sizeof ( Fluid* ), false );
  66         AddAttribute ( 0, "tag", sizeof ( bool ), false );
  67
  68         SPH_Setup ();
  69         Reset ( total );
  70 }
  71
  72 void FluidSystem::Reset ( int nmax )
  73 {
  74         ResetBuffer ( 0, nmax );
  75
  76         m_DT = 0.003; //  0.001;                        // .001 = for point grav
  77
  78         // Reset parameters
  79         m_Param [ MAX_FRAC ] = 1.0;
  80         m_Param [ POINT_GRAV ] = 0.0;
  81         m_Param [ PLANE_GRAV ] = 1.0;
  82
  83         m_Param [ BOUND_ZMIN_SLOPE ] = 0.0;
  84         m_Param [ FORCE_XMAX_SIN ] = 0.0;
  85         m_Param [ FORCE_XMIN_SIN ] = 0.0;
  86         m_Toggle [ WRAP_X ] = false;
  87         m_Toggle [ WALL_BARRIER ] = false;
  88         m_Toggle [ LEVY_BARRIER ] = false;
  89         m_Toggle [ DRAIN_BARRIER ] = false;
  90         m_Param [ SPH_INTSTIFF ] = 1.00;
  91         m_Param [ SPH_VISC ] = 0.2;
  92         m_Param [ SPH_INTSTIFF ] = 0.50;
  93         m_Param [ SPH_EXTSTIFF ] = 20000;
  94         m_Param [ SPH_SMOOTHRADIUS ] = 0.01;
  95
  96         m_Vec [ POINT_GRAV_POS ].Set ( 0, 0, 50 );
  97         m_Vec [ PLANE_GRAV_DIR ].Set ( 0, 0, -9.8 );
  98         m_Vec [ EMIT_POS ].Set ( 0, 0, 0 );
  99         m_Vec [ EMIT_RATE ].Set ( 0, 0, 0 );
 100         m_Vec [ EMIT_ANG ].Set ( 0, 90, 1.0 );
 101         m_Vec [ EMIT_DANG ].Set ( 0, 0, 0 );
 102 }
 103
 104 int FluidSystem::AddPoint ()
 105 {
 106         xref ndx;
 107         Fluid* f = (Fluid*) AddElem ( 0, ndx );
 108         f->sph_force.Set(0,0,0);
 109         f->vel.Set(0,0,0);
 110         f->vel_eval.Set(0,0,0);
 111         f->next = 0x0;
 112         f->pressure = 0;
 113         f->density = 0;
 114         return ndx;
 115 }
 116
 117 int FluidSystem::AddPointReuse ()
 118 {
 119         xref ndx;
 120         Fluid* f;
 121         if ( NumPoints() <= mBuf[0].max-2 )
 122                 f = (Fluid*) AddElem ( 0, ndx );
 123         else
 124                 f = (Fluid*) RandomElem ( 0, ndx );
 125
 126         f->sph_force.Set(0,0,0);
 127         f->vel.Set(0,0,0);
 128         f->vel_eval.Set(0,0,0);
 129         f->next = 0x0;
 130         f->pressure = 0;
 131         f->density = 0;
 132         return ndx;
 133 }
 134
 135 void FluidSystem::Run ()
 136 {
 137         bool bTiming = false;//true;
 138
 139         mint::Time start, stop;
 140
 141         float ss = m_Param [ SPH_PDIST ] / m_Param[ SPH_SIMSCALE ];             // simulation scale (not Schutzstaffel)
 142
 143         if ( m_Vec[EMIT_RATE].x > 0 && (++m_Frame) % (int) m_Vec[EMIT_RATE].x == 0 ) {
 144                 //m_Frame = 0;
 145                 Emit ( ss );
 146         }
 147
 148         #ifdef NOGRID
 149                 // Slow method - O(n^2)
 150                 SPH_ComputePressureSlow ();
 151                 SPH_ComputeForceSlow ();
 152         #else
 153
 154                 if ( m_Toggle[USE_CUDA] ) {
 155
 156                         #ifdef BUILD_CUDA
 157                                 // -- GPU --
 158                                 start.SetSystemTime ( ACC_NSEC );
 159                                 TransferToCUDA ( mBuf[0].data, (int*) &m_Grid[0], NumPoints() );
 160                                 if ( bTiming) { stop.SetSystemTime ( ACC_NSEC ); stop = stop - start; printf ( "TO: %s\n", stop.GetReadableTime().c_str() ); }
 161
 162                                 start.SetSystemTime ( ACC_NSEC );
 163                                 Grid_InsertParticlesCUDA ();
 164                                 if ( bTiming) { stop.SetSystemTime ( ACC_NSEC ); stop = stop - start; printf ( "INSERT (CUDA): %s\n", stop.GetReadableTime().c_str() ); }
 165
 166                                 start.SetSystemTime ( ACC_NSEC );
 167                                 SPH_ComputePressureCUDA ();
 168                                 if ( bTiming) { stop.SetSystemTime ( ACC_NSEC ); stop = stop - start; printf ( "PRESS (CUDA): %s\n", stop.GetReadableTime().c_str() ); }
 169
 170                                 start.SetSystemTime ( ACC_NSEC );
 171                                 SPH_ComputeForceCUDA ();
 172                                 if ( bTiming) { stop.SetSystemTime ( ACC_NSEC ); stop = stop - start; printf ( "FORCE (CUDA): %s\n", stop.GetReadableTime().c_str() ); }
 173
 174                                 //** CUDA integrator is incomplete..
 175                                 // Once integrator is done, we can remove TransferTo/From steps
 176                                 /*start.SetSystemTime ( ACC_NSEC );
 177                                 SPH_AdvanceCUDA( m_DT, m_DT/m_Param[SPH_SIMSCALE] );
 178                                 if ( bTiming) { stop.SetSystemTime ( ACC_NSEC ); stop = stop - start; printf ( "ADV (CUDA): %s\n", stop.GetReadableTime().c_str() ); }*/
 179
 180                                 start.SetSystemTime ( ACC_NSEC );
 181                                 TransferFromCUDA ( mBuf[0].data, (int*) &m_Grid[0], NumPoints() );
 182                                 if ( bTiming) { stop.SetSystemTime ( ACC_NSEC ); stop = stop - start; printf ( "FROM: %s\n", stop.GetReadableTime().c_str() ); }
 183
 184                                 // .. Do advance on CPU
 185                                 Advance();
 186
 187                         #endif
 188
 189                 } else {
 190                         // -- CPU only --
 191
 192                         start.SetSystemTime ( ACC_NSEC );
 193                         Grid_InsertParticles ();
 194                         if ( bTiming) { stop.SetSystemTime ( ACC_NSEC ); stop = stop - start; printf ( "INSERT: %s\n", stop.GetReadableTime().c_str() ); }
 195
 196                         start.SetSystemTime ( ACC_NSEC );
 197                         SPH_ComputePressureGrid ();
 198                         if ( bTiming) { stop.SetSystemTime ( ACC_NSEC ); stop = stop - start; printf ( "PRESS: %s\n", stop.GetReadableTime().c_str() ); }
 199
 200                         start.SetSystemTime ( ACC_NSEC );
 201                         SPH_ComputeForceGridNC ();
 202                         if ( bTiming) { stop.SetSystemTime ( ACC_NSEC ); stop = stop - start; printf ( "FORCE: %s\n", stop.GetReadableTime().c_str() ); }
 203
 204                         start.SetSystemTime ( ACC_NSEC );
 205                         Advance();
 206                         if ( bTiming) { stop.SetSystemTime ( ACC_NSEC ); stop = stop - start; printf ( "ADV: %s\n", stop.GetReadableTime().c_str() ); }
 207                 }
 208
 209         #endif
 210 }
 211
 212
 213
 214 void FluidSystem::SPH_DrawDomain ()
 215 {
 216         Vector3DF min, max;
 217         min = m_Vec[SPH_VOLMIN];
 218         max = m_Vec[SPH_VOLMAX];
 219         min.z += 0.5;
 220
 221         glColor3f ( 0.0, 0.0, 1.0 );
 222         glBegin ( GL_LINES );
 223         glVertex3f ( min.x, min.y, min.z );     glVertex3f ( max.x, min.y, min.z );
 224         glVertex3f ( min.x, max.y, min.z );     glVertex3f ( max.x, max.y, min.z );
 225         glVertex3f ( min.x, min.y, min.z );     glVertex3f ( min.x, max.y, min.z );
 226         glVertex3f ( max.x, min.y, min.z );     glVertex3f ( max.x, max.y, min.z );
 227         glEnd ();
 228 }
 229
 230 void FluidSystem::Advance ()
 231 {
 232         char *dat1, *dat1_end;
 233         Fluid* p;
 234         Vector3DF norm, z;
 235         Vector3DF dir, accel;
 236         Vector3DF vnext;
 237         Vector3DF min, max;
 238         double adj;
 239         float SL, SL2, ss, radius;
 240         float stiff, damp, speed, diff;
 241         SL = m_Param[SPH_LIMIT];
 242         SL2 = SL*SL;
 243
 244         stiff = m_Param[SPH_EXTSTIFF];
 245         damp = m_Param[SPH_EXTDAMP];
 246         radius = m_Param[SPH_PRADIUS];
 247         min = m_Vec[SPH_VOLMIN];
 248         max = m_Vec[SPH_VOLMAX];
 249         ss = m_Param[SPH_SIMSCALE];
 250
 251         dat1_end = mBuf[0].data + NumPoints()*mBuf[0].stride;
 252         for ( dat1 = mBuf[0].data; dat1 < dat1_end; dat1 += mBuf[0].stride ) {
 253                 p = (Fluid*) dat1;
 254
 255                 // Compute Acceleration
 256                 accel = p->sph_force;
 257                 accel *= m_Param[SPH_PMASS];
 258
 259                 // Velocity limiting
 260                 speed = accel.x*accel.x + accel.y*accel.y + accel.z*accel.z;
 261                 if ( speed > SL2 ) {
 262                         accel *= SL / sqrt(speed);
 263                 }
 264
 265                 // Boundary Conditions
 266
 267                 // Z-axis walls
 268                 diff = 2 * radius - ( p->pos.z - min.z - (p->pos.x - m_Vec[SPH_VOLMIN].x) * m_Param[BOUND_ZMIN_SLOPE] )*ss;
 269                 if (diff > EPSILON ) {
 270                         norm.Set ( -m_Param[BOUND_ZMIN_SLOPE], 0, 1.0 - m_Param[BOUND_ZMIN_SLOPE] );
 271                         adj = stiff * diff - damp * norm.Dot ( p->vel_eval );
 272                         accel.x += adj * norm.x; accel.y += adj * norm.y; accel.z += adj * norm.z;
 273                 }
 274
 275                 diff = 2 * radius - ( max.z - p->pos.z )*ss;
 276                 if (diff > EPSILON) {
 277                         norm.Set ( 0, 0, -1 );
 278                         adj = stiff * diff - damp * norm.Dot ( p->vel_eval );
 279                         accel.x += adj * norm.x; accel.y += adj * norm.y; accel.z += adj * norm.z;
 280                 }
 281
 282                 // X-axis walls
 283                 if ( !m_Toggle[WRAP_X] ) {
 284                         diff = 2 * radius - ( p->pos.x - min.x + (sin(m_Time*10.0)-1+(p->pos.y*0.025)*0.25) * m_Param[FORCE_XMIN_SIN] )*ss;
 285                         //diff = 2 * radius - ( p->pos.x - min.x + (sin(m_Time*10.0)-1) * m_Param[FORCE_XMIN_SIN] )*ss;
 286                         if (diff > EPSILON ) {
 287                                 norm.Set ( 1.0, 0, 0 );
 288                                 adj = (m_Param[ FORCE_XMIN_SIN ] + 1) * stiff * diff - damp * norm.Dot ( p->vel_eval ) ;
 289                                 accel.x += adj * norm.x; accel.y += adj * norm.y; accel.z += adj * norm.z;
 290                         }
 291
 292                         diff = 2 * radius - ( max.x - p->pos.x + (sin(m_Time*10.0)-1) * m_Param[FORCE_XMAX_SIN] )*ss;
 293                         if (diff > EPSILON) {
 294                                 norm.Set ( -1, 0, 0 );
 295                                 adj = (m_Param[ FORCE_XMAX_SIN ]+1) * stiff * diff - damp * norm.Dot ( p->vel_eval );
 296                                 accel.x += adj * norm.x; accel.y += adj * norm.y; accel.z += adj * norm.z;
 297                         }
 298                 }
 299
 300                 // Y-axis walls
 301                 diff = 2 * radius - ( p->pos.y - min.y )*ss;
 302                 if (diff > EPSILON) {
 303                         norm.Set ( 0, 1, 0 );
 304                         adj = stiff * diff - damp * norm.Dot ( p->vel_eval );
 305                         accel.x += adj * norm.x; accel.y += adj * norm.y; accel.z += adj * norm.z;
 306                 }
 307                 diff = 2 * radius - ( max.y - p->pos.y )*ss;
 308                 if (diff > EPSILON) {
 309                         norm.Set ( 0, -1, 0 );
 310                         adj = stiff * diff - damp * norm.Dot ( p->vel_eval );
 311                         accel.x += adj * norm.x; accel.y += adj * norm.y; accel.z += adj * norm.z;
 312                 }
 313
 314                 // Wall barrier
 315                 if ( m_Toggle[WALL_BARRIER] ) {
 316                         diff = 2 * radius - ( p->pos.x - 0 )*ss;
 317                         if (diff < 2*radius && diff > EPSILON && fabs(p->pos.y) < 3 && p->pos.z < 10) {
 318                                 norm.Set ( 1.0, 0, 0 );
 319                                 adj = 2*stiff * diff - damp * norm.Dot ( p->vel_eval ) ;
 320                                 accel.x += adj * norm.x; accel.y += adj * norm.y; accel.z += adj * norm.z;
 321                         }
 322                 }
 323
 324                 // Levy barrier
 325                 if ( m_Toggle[LEVY_BARRIER] ) {
 326                         diff = 2 * radius - ( p->pos.x - 0 )*ss;
 327                         if (diff < 2*radius && diff > EPSILON && fabs(p->pos.y) > 5 && p->pos.z < 10) {
 328                                 norm.Set ( 1.0, 0, 0 );
 329                                 adj = 2*stiff * diff - damp * norm.Dot ( p->vel_eval ) ;
 330                                 accel.x += adj * norm.x; accel.y += adj * norm.y; accel.z += adj * norm.z;
 331                         }
 332                 }
 333                 // Drain barrier
 334                 if ( m_Toggle[DRAIN_BARRIER] ) {
 335                         diff = 2 * radius - ( p->pos.z - min.z-15 )*ss;
 336                         if (diff < 2*radius && diff > EPSILON && (fabs(p->pos.x)>3 || fabs(p->pos.y)>3) ) {
 337                                 norm.Set ( 0, 0, 1);
 338                                 adj = stiff * diff - damp * norm.Dot ( p->vel_eval );
 339                                 accel.x += adj * norm.x; accel.y += adj * norm.y; accel.z += adj * norm.z;
 340                         }
 341                 }
 342
 343                 // Plane gravity
 344                 if ( m_Param[PLANE_GRAV] > 0)
 345                         accel += m_Vec[PLANE_GRAV_DIR];
 346
 347                 // Point gravity
 348                 if ( m_Param[POINT_GRAV] > 0 ) {
 349                         norm.x = ( p->pos.x - m_Vec[POINT_GRAV_POS].x );
 350                         norm.y = ( p->pos.y - m_Vec[POINT_GRAV_POS].y );
 351                         norm.z = ( p->pos.z - m_Vec[POINT_GRAV_POS].z );
 352                         norm.Normalize ();
 353                         norm *= m_Param[POINT_GRAV];
 354                         accel -= norm;
 355                 }
 356
 357                 // Leapfrog Integration ----------------------------
 358                 vnext = accel;
 359                 vnext *= m_DT;
 360                 vnext += p->vel;                                                // v(t+1/2) = v(t-1/2) + a(t) dt
 361                 p->vel_eval = p->vel;
 362                 p->vel_eval += vnext;
 363                 p->vel_eval *= 0.5;                                     // v(t+1) = [v(t-1/2) + v(t+1/2)] * 0.5         used to compute forces later
 364                 p->vel = vnext;
 365                 vnext *= m_DT/ss;
 366                 p->pos += vnext;                                                // p(t+1) = p(t) + v(t+1/2) dt
 367
 368                 if ( m_Param[CLR_MODE]==1.0 ) {
 369                         adj = fabs(vnext.x)+fabs(vnext.y)+fabs(vnext.z) / 7000.0;
 370                         adj = (adj > 1.0) ? 1.0 : adj;
 371                         p->clr = COLORA( 0, adj, adj, 1 );
 372                 }
 373                 if ( m_Param[CLR_MODE]==2.0 ) {
 374                         float v = 0.5 + ( p->pressure / 1500.0);
 375                         if ( v < 0.1 ) v = 0.1;
 376                         if ( v > 1.0 ) v = 1.0;
 377                         p->clr = COLORA ( v, 1-v, 0, 1 );
 378                 }
 379
 380
 381                 // Euler integration -------------------------------
 382                 /* accel += m_Gravity;
 383                 accel *= m_DT;
 384                 p->vel += accel;                                // v(t+1) = v(t) + a(t) dt
 385                 p->vel_eval += accel;
 386                 p->vel_eval *= m_DT/d;
 387                 p->pos += p->vel_eval;
 388                 p->vel_eval = p->vel;  */
 389
 390
 391                 if ( m_Toggle[WRAP_X] ) {
 392                         diff = p->pos.x - (m_Vec[SPH_VOLMIN].x + 2);                    // -- Simulates object in center of flow
 393                         if ( diff <= 0 ) {
 394                                 p->pos.x = (m_Vec[SPH_VOLMAX].x - 2) + diff*2;
 395                                 p->pos.z = 10;
 396                         }
 397                 }
 398         }
 399
 400         m_Time += m_DT;
 401 }
 402
 403 //------------------------------------------------------ SPH Setup
 404 //
 405 //  Range = +/- 10.0 * 0.006 (r) =         0.12                 m (= 120 mm = 4.7 inch)
 406 //  Container Volume (Vc) =                        0.001728             m^3
 407 //  Rest Density (D) =                          1000.0                  kg / m^3
 408 //  Particle Mass (Pm) =                           0.00020543   kg                                              (mass = vol * density)
 409 //  Number of Particles (N) =           4000.0
 410 //  Water Mass (M) =                               0.821                kg (= 821 grams)
 411 //  Water Volume (V) =                             0.000821     m^3 (= 3.4 cups, .21 gals)
 412 //  Smoothing Radius (R) =             0.02                     m (= 20 mm = ~3/4 inch)
 413 //  Particle Radius (Pr) =                         0.00366              m (= 4 mm  = ~1/8 inch)
 414 //  Particle Volume (Pv) =                         2.054e-7             m^3     (= .268 milliliters)
 415 //  Rest Distance (Pd) =                           0.0059               m
 416 //
 417 //  Given: D, Pm, N
 418 //    Pv = Pm / D                       0.00020543 kg / 1000 kg/m^3 = 2.054e-7 m^3
 419 //    Pv = 4/3*pi*Pr^3    cuberoot( 2.054e-7 m^3 * 3/(4pi) ) = 0.00366 m
 420 //     M = Pm * N                       0.00020543 kg * 4000.0 = 0.821 kg
 421 //     V =  M / D              0.821 kg / 1000 kg/m^3 = 0.000821 m^3
 422 //     V = Pv * N                        2.054e-7 m^3 * 4000 = 0.000821 m^3
 423 //    Pd = cuberoot(Pm/D)    cuberoot(0.00020543/1000) = 0.0059 m
 424 //
 425 // Ideal grid cell size (gs) = 2 * smoothing radius = 0.02*2 = 0.04
 426 // Ideal domain size = k*gs/d = k*0.02*2/0.005 = k*8 = {8, 16, 24, 32, 40, 48, ..}
 427 //    (k = number of cells, gs = cell size, d = simulation scale)
 428
 429 void FluidSystem::SPH_Setup ()
 430 {
 431         m_Param [ SPH_SIMSCALE ] =              0.004;                  // unit size
 432         m_Param [ SPH_VISC ] =                  0.2;                    // pascal-second (Pa.s) = 1 kg m^-1 s^-1  (see wikipedia page on viscosity)
 433         m_Param [ SPH_RESTDENSITY ] =   600.0;                  // kg / m^3
 434         m_Param [ SPH_PMASS ] =                 0.00020543;             // kg
 435         m_Param [ SPH_PRADIUS ] =               0.004;                  // m
 436         m_Param [ SPH_PDIST ] =                 0.0059;                 // m
 437         m_Param [ SPH_SMOOTHRADIUS ] =  0.01;                   // m
 438         m_Param [ SPH_INTSTIFF ] =              1.00;
 439         m_Param [ SPH_EXTSTIFF ] =              10000.0;
 440         m_Param [ SPH_EXTDAMP ] =               256.0;
 441         m_Param [ SPH_LIMIT ] =                 200.0;                  // m / s
 442
 443         m_Toggle [ SPH_GRID ] =         false;
 444         m_Toggle [ SPH_DEBUG ] =        false;
 445
 446         SPH_ComputeKernels ();
 447 }
 448
 449 void FluidSystem::SPH_ComputeKernels ()
 450 {
 451         m_Param [ SPH_PDIST ] = pow ( m_Param[SPH_PMASS] / m_Param[SPH_RESTDENSITY], 1/3.0 );
 452         m_R2 = m_Param [SPH_SMOOTHRADIUS] * m_Param[SPH_SMOOTHRADIUS];
 453         m_Poly6Kern = 315.0f / (64.0f * 3.141592 * pow( m_Param[SPH_SMOOTHRADIUS], 9) );        // Wpoly6 kernel (denominator part) - 2003 Muller, p.4
 454         m_SpikyKern = -45.0f / (3.141592 * pow( m_Param[SPH_SMOOTHRADIUS], 6) );                        // Laplacian of viscocity (denominator): PI h^6
 455         m_LapKern = 45.0f / (3.141592 * pow( m_Param[SPH_SMOOTHRADIUS], 6) );
 456 }
 457
 458 void FluidSystem::SPH_CreateExample ( int n, int nmax )
 459 {
 460         Vector3DF pos;
 461         Vector3DF min, max;
 462
 463         Reset ( nmax );
 464
 465         switch ( n ) {
 466         case 0:         // Wave pool
 467
 468                 //-- TEST CASE: 2x2x2 grid, 32 particles.  NOTE: Set PRADIUS to 0.0004 to reduce wall influence
 469                 //     grid 0:    3*3*2 = 18 particles
 470                 //     grid 1,2:  3*1*2 =  6 particles
 471                 //     grid 3:    1*1*2 =  2 particles
 472                 //     grid 4,5,6:    0 =  0 particles
 473                 /*m_Vec [ SPH_VOLMIN ].Set ( -2.5, -2.5, 0 );
 474                 m_Vec [ SPH_VOLMAX ].Set ( 2.5, 2.5, 5.0 );
 475                 m_Vec [ SPH_INITMIN ].Set ( -2.5, -2.5, 0 );
 476                 m_Vec [ SPH_INITMAX ].Set ( 2.5, 2.5, 1.6 );*/
 477
 478                 m_Vec [ SPH_VOLMIN ].Set ( -30, -30, 0 );
 479                 m_Vec [ SPH_VOLMAX ].Set ( 30, 30, 40 );
 480
 481                 //m_Vec [ SPH_INITMIN ].Set ( -5, -5, 10 );
 482                 //m_Vec [ SPH_INITMAX ].Set ( 5, 5, 20 );
 483
 484                 m_Vec [ SPH_INITMIN ].Set ( -20, -26, 10 );
 485                 m_Vec [ SPH_INITMAX ].Set ( 20, 26, 40 );
 486
 487                 m_Param [ FORCE_XMIN_SIN ] = 12.0;
 488                 m_Param [ BOUND_ZMIN_SLOPE ] = 0.05;
 489                 break;
 490         case 1:         // Dam break
 491                 m_Vec [ SPH_VOLMIN ].Set ( -30, -14, 0 );
 492                 m_Vec [ SPH_VOLMAX ].Set ( 30, 14, 60 );
 493                 m_Vec [ SPH_INITMIN ].Set ( 0, -13, 0 );
 494                 m_Vec [ SPH_INITMAX ].Set ( 29, 13, 30 );
 495                 m_Vec [ PLANE_GRAV_DIR ].Set ( 0.0, 0, -9.8 );
 496                 break;
 497         case 2:         // Dual-Wave pool
 498                 m_Vec [ SPH_VOLMIN ].Set ( -60, -5, 0 );
 499                 m_Vec [ SPH_VOLMAX ].Set ( 60, 5, 50 );
 500                 m_Vec [ SPH_INITMIN ].Set ( -46, -5, 0 );
 501                 m_Vec [ SPH_INITMAX ].Set ( 46, 5, 15 );
 502                 m_Param [ FORCE_XMIN_SIN ] = 8.0;
 503                 m_Param [ FORCE_XMAX_SIN ] = 8.0;
 504                 m_Vec [ PLANE_GRAV_DIR ].Set ( 0.0, 0, -9.8 );
 505                 break;
 506         case 3:         // Swirl Stream
 507                 m_Vec [ SPH_VOLMIN ].Set ( -30, -30, 0 );
 508                 m_Vec [ SPH_VOLMAX ].Set ( 30, 30, 50 );
 509                 m_Vec [ SPH_INITMIN ].Set ( -30, -30, 0 );
 510                 m_Vec [ SPH_INITMAX ].Set ( 30, 30, 40 );
 511                 m_Vec [ EMIT_POS ].Set ( -20, -20, 22 );
 512                 m_Vec [ EMIT_RATE ].Set ( 1, 4, 0 );
 513                 m_Vec [ EMIT_ANG ].Set ( 0, 120, 1.5 );
 514                 m_Vec [ EMIT_DANG ].Set ( 0, 0, 0 );
 515                 m_Vec [ PLANE_GRAV_DIR ].Set ( 0.0, 0, -9.8 );
 516                 break;
 517         case 4:         // Shockwave
 518                 m_Vec [ SPH_VOLMIN ].Set ( -60, -15, 0 );
 519                 m_Vec [ SPH_VOLMAX ].Set ( 60, 15, 50 );
 520                 m_Vec [ SPH_INITMIN ].Set ( -59, -14, 0 );
 521                 m_Vec [ SPH_INITMAX ].Set ( 59, 14, 30 );
 522                 m_Vec [ PLANE_GRAV_DIR ].Set ( 0.0, 0, -9.8 );
 523                 m_Toggle [ WALL_BARRIER ] = true;
 524                 m_Toggle [ WRAP_X ] = true;
 525                 break;
 526         case 5:         // Zero gravity
 527                 m_Vec [ SPH_VOLMIN ].Set ( -40, -40, 0 );
 528                 m_Vec [ SPH_VOLMAX ].Set ( 40, 40, 50 );
 529                 m_Vec [ SPH_INITMIN ].Set ( -20, -20, 20 );
 530                 m_Vec [ SPH_INITMAX ].Set ( 20, 20, 40 );
 531                 m_Vec [ EMIT_POS ].Set ( -20, 0, 40 );
 532                 m_Vec [ EMIT_RATE ].Set ( 2, 1, 0 );
 533                 m_Vec [ EMIT_ANG ].Set ( 0, 120, 0.25 );
 534                 m_Vec [ PLANE_GRAV_DIR ].Set ( 0, 0, 0 );
 535                 m_Param [ SPH_INTSTIFF ] = 0.20;
 536                 break;
 537         case 6:         // Point gravity
 538                 m_Vec [ SPH_VOLMIN ].Set ( -40, -40, 0 );
 539                 m_Vec [ SPH_VOLMAX ].Set ( 40, 40, 50 );
 540                 m_Vec [ SPH_INITMIN ].Set ( -20, -20, 20 );
 541                 m_Vec [ SPH_INITMAX ].Set ( 20, 20, 40 );
 542                 m_Param [ SPH_INTSTIFF ] = 0.50;
 543                 m_Vec [ EMIT_POS ].Set ( -20, 20, 25 );
 544                 m_Vec [ EMIT_RATE ].Set ( 1, 4, 0 );
 545                 m_Vec [ EMIT_ANG ].Set ( -20, 100, 2.0 );
 546                 m_Vec [ POINT_GRAV_POS ].Set ( 0, 0, 25 );
 547                 m_Vec [ PLANE_GRAV_DIR ].Set ( 0, 0, 0 );
 548                 m_Param [ POINT_GRAV ] = 3.5;
 549                 break;
 550         case 7:         // Levy break
 551                 m_Vec [ SPH_VOLMIN ].Set ( -40, -40, 0 );
 552                 m_Vec [ SPH_VOLMAX ].Set ( 40, 40, 50 );
 553                 m_Vec [ SPH_INITMIN ].Set ( 10, -40, 0 );
 554                 m_Vec [ SPH_INITMAX ].Set ( 40, 40, 50 );
 555                 m_Vec [ EMIT_POS ].Set ( 34, 27, 16.6 );
 556                 m_Vec [ EMIT_RATE ].Set ( 2, 9, 0 );
 557                 m_Vec [ EMIT_ANG ].Set ( 118, 200, 1.0 );
 558                 m_Toggle [ LEVY_BARRIER ] = true;
 559                 m_Param [ BOUND_ZMIN_SLOPE ] = 0.1;
 560                 m_Vec [ PLANE_GRAV_DIR ].Set ( 0.0, 0, -9.8 );
 561                 break;
 562         case 8:         // Drain
 563                 m_Vec [ SPH_VOLMIN ].Set ( -20, -20, 0 );
 564                 m_Vec [ SPH_VOLMAX ].Set ( 20, 20, 50 );
 565                 m_Vec [ SPH_INITMIN ].Set ( -15, -20, 20 );
 566                 m_Vec [ SPH_INITMAX ].Set ( 20, 20, 50 );
 567                 m_Vec [ EMIT_POS ].Set ( -16, -16, 30 );
 568                 m_Vec [ EMIT_RATE ].Set ( 1, 4, 0 );
 569                 m_Vec [ EMIT_ANG ].Set ( -20, 140, 1.8 );
 570                 m_Toggle [ DRAIN_BARRIER ] = true;
 571                 m_Vec [ PLANE_GRAV_DIR ].Set ( 0.0, 0, -9.8 );
 572                 break;
 573         case 9:                 // Tumbler
 574                 m_Vec [ SPH_VOLMIN ].Set ( -30, -30, 0 );
 575                 m_Vec [ SPH_VOLMAX ].Set ( 30, 30, 50 );
 576                 m_Vec [ SPH_INITMIN ].Set ( 24, -29, 20 );
 577                 m_Vec [ SPH_INITMAX ].Set ( 29, 29, 40 );
 578                 m_Param [ SPH_VISC ] = 0.1;
 579                 m_Param [ SPH_INTSTIFF ] = 0.50;
 580                 m_Param [ SPH_EXTSTIFF ] = 8000;
 581                 //m_Param [ SPH_SMOOTHRADIUS ] = 0.01;
 582                 m_Param [ BOUND_ZMIN_SLOPE ] = 0.4;
 583                 m_Param [ FORCE_XMIN_SIN ] = 12.00;
 584                 m_Vec [ PLANE_GRAV_DIR ].Set ( 0.0, 0, -9.8 );
 585                 break;
 586         case 10:                // Large sim
 587                 m_Vec [ SPH_VOLMIN ].Set ( -35, -35, 0 );
 588                 m_Vec [ SPH_VOLMAX ].Set ( 35, 35, 60 );
 589                 m_Vec [ SPH_INITMIN ].Set ( -5, -35, 0 );
 590                 m_Vec [ SPH_INITMAX ].Set ( 30, 0, 60 );
 591                 m_Vec [ PLANE_GRAV_DIR ].Set ( 0.0, 0, -9.8 );
 592                 break;
 593         }
 594
 595         SPH_ComputeKernels ();
 596
 597         m_Param [ SPH_SIMSIZE ] = m_Param [ SPH_SIMSCALE ] * (m_Vec[SPH_VOLMAX].z - m_Vec[SPH_VOLMIN].z);
 598         m_Param [ SPH_PDIST ] = pow ( m_Param[SPH_PMASS] / m_Param[SPH_RESTDENSITY], 1/3.0 );
 599
 600         float ss = m_Param [ SPH_PDIST ]*0.87 / m_Param[ SPH_SIMSCALE ];
 601         printf ( "Spacing: %f\n", ss);
 602         AddVolume ( m_Vec[SPH_INITMIN], m_Vec[SPH_INITMAX], ss );       // Create the particles
 603
 604         float cell_size = m_Param[SPH_SMOOTHRADIUS]*2.0;                        // Grid cell size (2r)
 605         Grid_Setup ( m_Vec[SPH_VOLMIN], m_Vec[SPH_VOLMAX], m_Param[SPH_SIMSCALE], cell_size, 1.0 );                                                                                             // Setup grid
 606         Grid_InsertParticles ();                                                                        // Insert particles
 607
 608         Vector3DF vmin, vmax;
 609         vmin =  m_Vec[SPH_VOLMIN];
 610         vmin -= Vector3DF(2,2,2);
 611         vmax =  m_Vec[SPH_VOLMAX];
 612         vmax += Vector3DF(2,2,-2);
 613
 614         #ifdef BUILD_CUDA
 615                 FluidClearCUDA ();
 616                 Sleep ( 500 );
 617
 618                 FluidSetupCUDA ( NumPoints(), sizeof(Fluid), *(float3*)& m_GridMin, *(float3*)& m_GridMax, *(float3*)& m_GridRes, *(float3*)& m_GridSize, (int) m_Vec[EMIT_RATE].x );
 619
 620                 Sleep ( 500 );
 621
 622                 FluidParamCUDA ( m_Param[SPH_SIMSCALE], m_Param[SPH_SMOOTHRADIUS], m_Param[SPH_PMASS], m_Param[SPH_RESTDENSITY], m_Param[SPH_INTSTIFF], m_Param[SPH_VISC] );
 623         #endif
 624
 625 }
 626
 627 // Compute Pressures - Very slow yet simple. O(n^2)
 628 void FluidSystem::SPH_ComputePressureSlow ()
 629 {
 630         char *dat1, *dat1_end;
 631         char *dat2, *dat2_end;
 632         Fluid *p, *q;
 633         int cnt = 0;
 634         double dx, dy, dz, sum, dsq, c;
 635         double d, d2, mR, mR2;
 636         d = m_Param[SPH_SIMSCALE];
 637         d2 = d*d;
 638         mR = m_Param[SPH_SMOOTHRADIUS];
 639         mR2 = mR*mR;
 640
 641         dat1_end = mBuf[0].data + NumPoints()*mBuf[0].stride;
 642         for ( dat1 = mBuf[0].data; dat1 < dat1_end; dat1 += mBuf[0].stride ) {
 643                 p = (Fluid*) dat1;
 644
 645                 sum = 0.0;
 646                 cnt = 0;
 647
 648                 dat2_end = mBuf[0].data + NumPoints()*mBuf[0].stride;
 649                 for ( dat2 = mBuf[0].data; dat2 < dat2_end; dat2 += mBuf[0].stride ) {
 650                         q = (Fluid*) dat2;
 651
 652                         if ( p==q ) continue;
 653                         dx = ( p->pos.x - q->pos.x)*d;          // dist in cm
 654                         dy = ( p->pos.y - q->pos.y)*d;
 655                         dz = ( p->pos.z - q->pos.z)*d;
 656                         dsq = (dx*dx + dy*dy + dz*dz);
 657                         if ( mR2 > dsq ) {
 658                                 c =  m_R2 - dsq;
 659                                 sum += c * c * c;
 660                                 cnt++;
 661                                 //if ( p == m_CurrP ) q->tag = true;
 662                         }
 663                 }
 664                 p->density = sum * m_Param[SPH_PMASS] * m_Poly6Kern ;
 665                 p->pressure = ( p->density - m_Param[SPH_RESTDENSITY] ) * m_Param[SPH_INTSTIFF];
 666                 p->density = 1.0f / p->density;
 667         }
 668 }
 669
 670 // Compute Pressures - Using spatial grid, and also create neighbor table
 671 void FluidSystem::SPH_ComputePressureGrid ()
 672 {
 673         char *dat1, *dat1_end;
 674         Fluid* p;
 675         Fluid* pcurr;
 676         int pndx;
 677         int i, cnt = 0;
 678         float dx, dy, dz, sum, dsq, c;
 679         float d, d2, mR, mR2;
 680         float radius = m_Param[SPH_SMOOTHRADIUS] / m_Param[SPH_SIMSCALE];
 681         d = m_Param[SPH_SIMSCALE];
 682         d2 = d*d;
 683         mR = m_Param[SPH_SMOOTHRADIUS];
 684         mR2 = mR*mR;
 685
 686         dat1_end = mBuf[0].data + NumPoints()*mBuf[0].stride;
 687         i = 0;
 688         for ( dat1 = mBuf[0].data; dat1 < dat1_end; dat1 += mBuf[0].stride, i++ ) {
 689                 p = (Fluid*) dat1;
 690
 691                 sum = 0.0;
 692                 m_NC[i] = 0;
 693
 694                 Grid_FindCells ( p->pos, radius );
 695                 for (int cell=0; cell < 8; cell++) {
 696                         if ( m_GridCell[cell] != -1 ) {
 697                                 pndx = m_Grid [ m_GridCell[cell] ];
 698                                 while ( pndx != -1 ) {
 699                                         pcurr = (Fluid*) (mBuf[0].data + pndx*mBuf[0].stride);
 700                                         if ( pcurr == p ) {pndx = pcurr->next; continue; }
 701                                         dx = ( p->pos.x - pcurr->pos.x)*d;              // dist in cm
 702                                         dy = ( p->pos.y - pcurr->pos.y)*d;
 703                                         dz = ( p->pos.z - pcurr->pos.z)*d;
 704                                         dsq = (dx*dx + dy*dy + dz*dz);
 705                                         if ( mR2 > dsq ) {
 706                                                 c =  m_R2 - dsq;
 707                                                 sum += c * c * c;
 708                                                 if ( m_NC[i] < MAX_NEIGHBOR ) {
 709                                                         m_Neighbor[i][ m_NC[i] ] = pndx;
 710                                                         m_NDist[i][ m_NC[i] ] = sqrt(dsq);
 711                                                         m_NC[i]++;
 712                                                 }
 713                                         }
 714                                         pndx = pcurr->next;
 715                                 }
 716                         }
 717                         m_GridCell[cell] = -1;
 718                 }
 719                 p->density = sum * m_Param[SPH_PMASS] * m_Poly6Kern ;
 720                 p->pressure = ( p->density - m_Param[SPH_RESTDENSITY] ) * m_Param[SPH_INTSTIFF];
 721                 p->density = 1.0f / p->density;
 722         }
 723 }
 724
 725 // Compute Forces - Very slow, but simple. O(n^2)
 726 void FluidSystem::SPH_ComputeForceSlow ()
 727 {
 728         char *dat1, *dat1_end;
 729         char *dat2, *dat2_end;
 730         Fluid *p, *q;
 731         Vector3DF force, fcurr;
 732         register double pterm, vterm, dterm;
 733         double c, r, d, sum, dsq;
 734         double dx, dy, dz;
 735         double mR, mR2, visc;
 736
 737         d = m_Param[SPH_SIMSCALE];
 738         mR = m_Param[SPH_SMOOTHRADIUS];
 739         mR2 = (mR*mR);
 740         visc = m_Param[SPH_VISC];
 741         vterm = m_LapKern * visc;
 742
 743         dat1_end = mBuf[0].data + NumPoints()*mBuf[0].stride;
 744         for ( dat1 = mBuf[0].data; dat1 < dat1_end; dat1 += mBuf[0].stride ) {
 745                 p = (Fluid*) dat1;
 746
 747                 sum = 0.0;
 748                 force.Set ( 0, 0, 0 );
 749
 750                 dat2_end = mBuf[0].data + NumPoints()*mBuf[0].stride;
 751                 for ( dat2 = mBuf[0].data; dat2 < dat2_end; dat2 += mBuf[0].stride ) {
 752                         q = (Fluid*) dat2;
 753
 754                         if ( p == q ) continue;
 755                         dx = ( p->pos.x - q->pos.x )*d;                 // dist in cm
 756                         dy = ( p->pos.y - q->pos.y )*d;
 757                         dz = ( p->pos.z - q->pos.z )*d;
 758                         dsq = (dx*dx + dy*dy + dz*dz);
 759                         if ( mR2 > dsq ) {
 760                                 r = sqrt ( dsq );
 761                                 c = (mR - r);
 762                                 pterm = -0.5f * c * m_SpikyKern * ( p->pressure + q->pressure) / r;
 763                                 dterm = c * p->density * q->density;
 764                                 force.x += ( pterm * dx + vterm * (q->vel_eval.x - p->vel_eval.x) ) * dterm;
 765                                 force.y += ( pterm * dy + vterm * (q->vel_eval.y - p->vel_eval.y) ) * dterm;
 766                                 force.z += ( pterm * dz + vterm * (q->vel_eval.z - p->vel_eval.z) ) * dterm;
 767                         }
 768                 }
 769                 p->sph_force = force;
 770         }
 771 }
 772
 773 // Compute Forces - Using spatial grid. Faster.
 774 void FluidSystem::SPH_ComputeForceGrid ()
 775 {
 776         char *dat1, *dat1_end;
 777         Fluid *p;
 778         Fluid *pcurr;
 779         int pndx;
 780         Vector3DF force, fcurr;
 781         register double pterm, vterm, dterm;
 782         double c, d, dsq, r;
 783         double dx, dy, dz;
 784         double mR, mR2, visc;
 785         float radius = m_Param[SPH_SMOOTHRADIUS] / m_Param[SPH_SIMSCALE];
 786
 787         d = m_Param[SPH_SIMSCALE];
 788         mR = m_Param[SPH_SMOOTHRADIUS];
 789         mR2 = (mR*mR);
 790         visc = m_Param[SPH_VISC];
 791
 792         dat1_end = mBuf[0].data + NumPoints()*mBuf[0].stride;
 793         for ( dat1 = mBuf[0].data; dat1 < dat1_end; dat1 += mBuf[0].stride ) {
 794                 p = (Fluid*) dat1;
 795
 796                 force.Set ( 0, 0, 0 );
 797
 798                 Grid_FindCells ( p->pos, radius );
 799                 for (int cell=0; cell < 8; cell++) {
 800                         if ( m_GridCell[cell] != -1 ) {
 801                                 pndx = m_Grid [ m_GridCell[cell] ];
 802                                 while ( pndx != -1 ) {
 803                                         pcurr = (Fluid*) (mBuf[0].data + pndx*mBuf[0].stride);
 804                                         if ( pcurr == p ) {pndx = pcurr->next; continue; }
 805
 806                                         dx = ( p->pos.x - pcurr->pos.x)*d;              // dist in cm
 807                                         dy = ( p->pos.y - pcurr->pos.y)*d;
 808                                         dz = ( p->pos.z - pcurr->pos.z)*d;
 809                                         dsq = (dx*dx + dy*dy + dz*dz);
 810                                         if ( mR2 > dsq ) {
 811                                                 r = sqrt ( dsq );
 812                                                 c = (mR - r);
 813                                                 pterm = -0.5f * c * m_SpikyKern * ( p->pressure + pcurr->pressure) / r;
 814                                                 dterm = c * p->density * pcurr->density;
 815                                                 vterm = m_LapKern * visc;
 816                                                 force.x += ( pterm * dx + vterm * (pcurr->vel_eval.x - p->vel_eval.x) ) * dterm;
 817                                                 force.y += ( pterm * dy + vterm * (pcurr->vel_eval.y - p->vel_eval.y) ) * dterm;
 818                                                 force.z += ( pterm * dz + vterm * (pcurr->vel_eval.z - p->vel_eval.z) ) * dterm;
 819                                         }
 820                                         pndx = pcurr->next;
 821                                 }
 822                         }
 823                 }
 824                 p->sph_force = force;
 825         }
 826 }
 827
 828 // Compute Forces - Using spatial grid with saved neighbor table. Fastest.
 829 void FluidSystem::SPH_ComputeForceGridNC ()
 830 {
 831         char *dat1, *dat1_end;
 832         Fluid *p;
 833         Fluid *pcurr;
 834         Vector3DF force, fcurr;
 835         register float pterm, vterm, dterm;
 836         int i;
 837         float c, d;
 838         float dx, dy, dz;
 839         float mR, mR2, visc;
 840
 841         d = m_Param[SPH_SIMSCALE];
 842         mR = m_Param[SPH_SMOOTHRADIUS];
 843         mR2 = (mR*mR);
 844         visc = m_Param[SPH_VISC];
 845
 846         dat1_end = mBuf[0].data + NumPoints()*mBuf[0].stride;
 847         i = 0;
 848
 849         for ( dat1 = mBuf[0].data; dat1 < dat1_end; dat1 += mBuf[0].stride, i++ ) {
 850                 p = (Fluid*) dat1;
 851
 852                 force.Set ( 0, 0, 0 );
 853                 for (int j=0; j < m_NC[i]; j++ ) {
 854                         pcurr = (Fluid*) (mBuf[0].data + m_Neighbor[i][j]*mBuf[0].stride);
 855                         dx = ( p->pos.x - pcurr->pos.x)*d;              // dist in cm
 856                         dy = ( p->pos.y - pcurr->pos.y)*d;
 857                         dz = ( p->pos.z - pcurr->pos.z)*d;
 858                         c = ( mR - m_NDist[i][j] );
 859                         pterm = -0.5f * c * m_SpikyKern * ( p->pressure + pcurr->pressure) / m_NDist[i][j];
 860                         dterm = c * p->density * pcurr->density;
 861                         vterm = m_LapKern * visc;
 862                         force.x += ( pterm * dx + vterm * (pcurr->vel_eval.x - p->vel_eval.x) ) * dterm;
 863                         force.y += ( pterm * dy + vterm * (pcurr->vel_eval.y - p->vel_eval.y) ) * dterm;
 864                         force.z += ( pterm * dz + vterm * (pcurr->vel_eval.z - p->vel_eval.z) ) * dterm;
 865                 }
 866                 p->sph_force = force;
 867         }
 868 }
 869