| //Code written by Richard O. Lee and Christian Bienia |
| //Modified by Christian Fensch |
| |
| #include <cstdlib> |
| #include <cstring> |
| |
| #include <iostream> |
| #include <fstream> |
| #include <math.h> |
| #include <assert.h> |
| |
| #include "fluid.hpp" |
| #include "cellpool.hpp" |
| |
| #ifdef ENABLE_VISUALIZATION |
| #include "fluidview.hpp" |
| #endif |
| |
| #ifdef ENABLE_PARSEC_HOOKS |
| #include <hooks.h> |
| #endif |
| |
| //Uncomment to add code to check that Courant–Friedrichs–Lewy condition is satisfied at runtime |
| //#define ENABLE_CFL_CHECK |
| |
| //Define ENABLE_STATISTICS to collect additional information about the particles at runtime. |
| //#define ENABLE_STATISTICS |
| |
| //////////////////////////////////////////////////////////////////////////////// |
| |
| cellpool pool; |
| |
| fptype restParticlesPerMeter, h, hSq; |
| fptype densityCoeff, pressureCoeff, viscosityCoeff; |
| |
| int nx, ny, nz; // number of grid cells in each dimension |
| Vec3 delta; // cell dimensions |
| int numParticles = 0; |
| int numCells = 0; |
| Cell *cells = NULL; |
| Cell *cells2 = NULL; |
| int *cnumPars = 0; |
| int *cnumPars2 = 0; |
| Cell **last_cells = NULL; //helper array with pointers to last cell structure of "cells" array lists |
| #ifdef ENABLE_VISUALIZATION |
| Vec3 vMax(0.0,0.0,0.0); |
| Vec3 vMin(0.0,0.0,0.0); |
| #endif |
| //////////////////////////////////////////////////////////////////////////////// |
| |
| void InitSim(char const *fileName) |
| { |
| std::cout << "Loading file \"" << fileName << "\"..." << std::endl; |
| std::ifstream file(fileName, std::ios::binary); |
| if(!file) { |
| std::cerr << "Error opening file. Aborting." << std::endl; |
| exit(1); |
| } |
| |
| //Always use single precision float variables b/c file format uses single precision |
| float restParticlesPerMeter_le; |
| int numParticles_le; |
| file.read((char *)&restParticlesPerMeter_le, FILE_SIZE_FLOAT); |
| file.read((char *)&numParticles_le, FILE_SIZE_INT); |
| if(!isLittleEndian()) { |
| restParticlesPerMeter = bswap_float(restParticlesPerMeter_le); |
| numParticles = bswap_int32(numParticles_le); |
| } else { |
| restParticlesPerMeter = restParticlesPerMeter_le; |
| numParticles = numParticles_le; |
| } |
| cellpool_init(&pool, numParticles); |
| |
| h = kernelRadiusMultiplier / restParticlesPerMeter; |
| hSq = h*h; |
| |
| #ifndef ENABLE_DOUBLE_PRECISION |
| fptype coeff1 = 315.0 / (64.0*pi*powf(h,9.0)); |
| fptype coeff2 = 15.0 / (pi*powf(h,6.0)); |
| fptype coeff3 = 45.0 / (pi*powf(h,6.0)); |
| #else |
| fptype coeff1 = 315.0 / (64.0*pi*pow(h,9.0)); |
| fptype coeff2 = 15.0 / (pi*pow(h,6.0)); |
| fptype coeff3 = 45.0 / (pi*pow(h,6.0)); |
| #endif //ENABLE_DOUBLE_PRECISION |
| fptype particleMass = 0.5*doubleRestDensity / (restParticlesPerMeter*restParticlesPerMeter*restParticlesPerMeter); |
| densityCoeff = particleMass * coeff1; |
| pressureCoeff = 3.0*coeff2 * 0.50*stiffnessPressure * particleMass; |
| viscosityCoeff = viscosity * coeff3 * particleMass; |
| |
| Vec3 range = domainMax - domainMin; |
| nx = (int)(range.x / h); |
| ny = (int)(range.y / h); |
| nz = (int)(range.z / h); |
| assert(nx >= 1 && ny >= 1 && nz >= 1); |
| numCells = nx*ny*nz; |
| std::cout << "Number of cells: " << numCells << std::endl; |
| delta.x = range.x / nx; |
| delta.y = range.y / ny; |
| delta.z = range.z / nz; |
| assert(delta.x >= h && delta.y >= h && delta.z >= h); |
| |
| //make sure Cell structure is multiple of estiamted cache line size |
| assert(sizeof(Cell) % CACHELINE_SIZE == 0); |
| //make sure helper Cell structure is in sync with real Cell structure |
| #pragma warning( disable : 1684) // warning #1684: conversion from pointer to same-sized integral type (potential portability problem) |
| assert(offsetof(struct Cell_aux, padding) == offsetof(struct Cell, padding)); |
| |
| #if defined(WIN32) |
| cells = (struct Cell*)_aligned_malloc(sizeof(struct Cell) * numCells, CACHELINE_SIZE); |
| cells2 = (struct Cell*)_aligned_malloc(sizeof(struct Cell) * numCells, CACHELINE_SIZE); |
| cnumPars = (int*)_aligned_malloc(sizeof(int) * numCells, CACHELINE_SIZE); |
| cnumPars2 = (int*)_aligned_malloc(sizeof(int) * numCells, CACHELINE_SIZE); |
| last_cells = (struct Cell **)_aligned_malloc(sizeof(struct Cell *) * numCells, CACHELINE_SIZE); |
| assert((cells!=NULL) && (cells2!=NULL) && (cnumPars!=NULL) && (cnumPars2!=NULL) && (last_cells!=NULL)); |
| #else |
| int rv0 = posix_memalign((void **)(&cells), CACHELINE_SIZE, sizeof(struct Cell) * numCells); |
| int rv1 = posix_memalign((void **)(&cells2), CACHELINE_SIZE, sizeof(struct Cell) * numCells); |
| int rv2 = posix_memalign((void **)(&cnumPars), CACHELINE_SIZE, sizeof(int) * numCells); |
| int rv3 = posix_memalign((void **)(&cnumPars2), CACHELINE_SIZE, sizeof(int) * numCells); |
| int rv4 = posix_memalign((void **)(&last_cells), CACHELINE_SIZE, sizeof(struct Cell *) * numCells); |
| assert((rv0==0) && (rv1==0) && (rv2==0) && (rv3==0) && (rv4==0)); |
| #endif |
| |
| // because cells and cells2 are not allocated via new |
| // we construct them here |
| for(int i=0; i<numCells; ++i) |
| { |
| new (&cells[i]) Cell; |
| new (&cells2[i]) Cell; |
| } |
| |
| memset(cnumPars, 0, numCells*sizeof(int)); |
| |
| //Always use single precision float variables b/c file format uses single precision |
| float px, py, pz, hvx, hvy, hvz, vx, vy, vz; |
| for(int i = 0; i < numParticles; ++i) |
| { |
| file.read((char *)&px, FILE_SIZE_FLOAT); |
| file.read((char *)&py, FILE_SIZE_FLOAT); |
| file.read((char *)&pz, FILE_SIZE_FLOAT); |
| file.read((char *)&hvx, FILE_SIZE_FLOAT); |
| file.read((char *)&hvy, FILE_SIZE_FLOAT); |
| file.read((char *)&hvz, FILE_SIZE_FLOAT); |
| file.read((char *)&vx, FILE_SIZE_FLOAT); |
| file.read((char *)&vy, FILE_SIZE_FLOAT); |
| file.read((char *)&vz, FILE_SIZE_FLOAT); |
| if(!isLittleEndian()) { |
| px = bswap_float(px); |
| py = bswap_float(py); |
| pz = bswap_float(pz); |
| hvx = bswap_float(hvx); |
| hvy = bswap_float(hvy); |
| hvz = bswap_float(hvz); |
| vx = bswap_float(vx); |
| vy = bswap_float(vy); |
| vz = bswap_float(vz); |
| } |
| int ci = (int)(((fptype)px - domainMin.x) / delta.x); |
| int cj = (int)(((fptype)py - domainMin.y) / delta.y); |
| int ck = (int)(((fptype)pz - domainMin.z) / delta.z); |
| |
| if(ci < 0) ci = 0; else if(ci >= nx) ci = nx-1; |
| if(cj < 0) cj = 0; else if(cj >= ny) cj = ny-1; |
| if(ck < 0) ck = 0; else if(ck >= nz) ck = nz-1; |
| |
| int index = (ck*ny + cj)*nx + ci; |
| Cell *cell = &cells[index]; |
| |
| //go to last cell structure in list |
| int np = cnumPars[index]; |
| while(np > PARTICLES_PER_CELL) { |
| cell = cell->next; |
| np = np - PARTICLES_PER_CELL; |
| } |
| //add another cell structure if everything full |
| if( (np % PARTICLES_PER_CELL == 0) && (cnumPars[index] != 0) ) { |
| cell->next = cellpool_getcell(&pool); |
| cell = cell->next; |
| np = np - PARTICLES_PER_CELL; // np = 0; |
| } |
| |
| //add particle to cell |
| cell->p[np].x = px; |
| cell->p[np].y = py; |
| cell->p[np].z = pz; |
| cell->hv[np].x = hvx; |
| cell->hv[np].y = hvy; |
| cell->hv[np].z = hvz; |
| cell->v[np].x = vx; |
| cell->v[np].y = vy; |
| cell->v[np].z = vz; |
| #ifdef ENABLE_VISUALIZATION |
| vMin.x = std::min(vMin.x, cell->v[np].x); |
| vMax.x = std::max(vMax.x, cell->v[np].x); |
| vMin.y = std::min(vMin.y, cell->v[np].y); |
| vMax.y = std::max(vMax.y, cell->v[np].y); |
| vMin.z = std::min(vMin.z, cell->v[np].z); |
| vMax.z = std::max(vMax.z, cell->v[np].z); |
| #endif |
| ++cnumPars[index]; |
| } |
| |
| std::cout << "Number of particles: " << numParticles << std::endl; |
| } |
| |
| //////////////////////////////////////////////////////////////////////////////// |
| |
| void SaveFile(char const *fileName) |
| { |
| std::cout << "Saving file \"" << fileName << "\"..." << std::endl; |
| |
| std::ofstream file(fileName, std::ios::binary); |
| assert(file); |
| |
| //Always use single precision float variables b/c file format uses single precision |
| if(!isLittleEndian()) { |
| float restParticlesPerMeter_le; |
| int numParticles_le; |
| |
| restParticlesPerMeter_le = bswap_float((float)restParticlesPerMeter); |
| numParticles_le = bswap_int32(numParticles); |
| file.write((char *)&restParticlesPerMeter_le, FILE_SIZE_FLOAT); |
| file.write((char *)&numParticles_le, FILE_SIZE_INT); |
| } else { |
| file.write((char *)&restParticlesPerMeter, FILE_SIZE_FLOAT); |
| file.write((char *)&numParticles, FILE_SIZE_INT); |
| } |
| |
| int count = 0; |
| for(int i = 0; i < numCells; ++i) |
| { |
| Cell *cell = &cells[i]; |
| int np = cnumPars[i]; |
| for(int j = 0; j < np; ++j) |
| { |
| //Always use single precision float variables b/c file format uses single precision |
| float px, py, pz, hvx, hvy, hvz, vx,vy, vz; |
| if(!isLittleEndian()) { |
| px = bswap_float((float)(cell->p[j % PARTICLES_PER_CELL].x)); |
| py = bswap_float((float)(cell->p[j % PARTICLES_PER_CELL].y)); |
| pz = bswap_float((float)(cell->p[j % PARTICLES_PER_CELL].z)); |
| hvx = bswap_float((float)(cell->hv[j % PARTICLES_PER_CELL].x)); |
| hvy = bswap_float((float)(cell->hv[j % PARTICLES_PER_CELL].y)); |
| hvz = bswap_float((float)(cell->hv[j % PARTICLES_PER_CELL].z)); |
| vx = bswap_float((float)(cell->v[j % PARTICLES_PER_CELL].x)); |
| vy = bswap_float((float)(cell->v[j % PARTICLES_PER_CELL].y)); |
| vz = bswap_float((float)(cell->v[j % PARTICLES_PER_CELL].z)); |
| } else { |
| px = (float)(cell->p[j % PARTICLES_PER_CELL].x); |
| py = (float)(cell->p[j % PARTICLES_PER_CELL].y); |
| pz = (float)(cell->p[j % PARTICLES_PER_CELL].z); |
| hvx = (float)(cell->hv[j % PARTICLES_PER_CELL].x); |
| hvy = (float)(cell->hv[j % PARTICLES_PER_CELL].y); |
| hvz = (float)(cell->hv[j % PARTICLES_PER_CELL].z); |
| vx = (float)(cell->v[j % PARTICLES_PER_CELL].x); |
| vy = (float)(cell->v[j % PARTICLES_PER_CELL].y); |
| vz = (float)(cell->v[j % PARTICLES_PER_CELL].z); |
| } |
| file.write((char *)&px, FILE_SIZE_FLOAT); |
| file.write((char *)&py, FILE_SIZE_FLOAT); |
| file.write((char *)&pz, FILE_SIZE_FLOAT); |
| file.write((char *)&hvx, FILE_SIZE_FLOAT); |
| file.write((char *)&hvy, FILE_SIZE_FLOAT); |
| file.write((char *)&hvz, FILE_SIZE_FLOAT); |
| file.write((char *)&vx, FILE_SIZE_FLOAT); |
| file.write((char *)&vy, FILE_SIZE_FLOAT); |
| file.write((char *)&vz, FILE_SIZE_FLOAT); |
| ++count; |
| |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) { |
| cell = cell->next; |
| } |
| } |
| } |
| assert(count == numParticles); |
| } |
| |
| //////////////////////////////////////////////////////////////////////////////// |
| |
| void CleanUpSim() |
| { |
| // first return extended cells to cell pools |
| for(int i=0; i< numCells; ++i) |
| { |
| Cell& cell = cells[i]; |
| while(cell.next) |
| { |
| Cell *temp = cell.next; |
| cell.next = temp->next; |
| cellpool_returncell(&pool, temp); |
| } |
| } |
| // now return cell pools |
| cellpool_destroy(&pool); |
| |
| #if defined(WIN32) |
| _aligned_free(cells); |
| _aligned_free(cells2); |
| _aligned_free(cnumPars); |
| _aligned_free(cnumPars2); |
| _aligned_free(last_cells); |
| #else |
| free(cells); |
| free(cells2); |
| free(cnumPars); |
| free(cnumPars2); |
| free(last_cells); |
| #endif |
| } |
| |
| //////////////////////////////////////////////////////////////////////////////// |
| |
| void RebuildGrid() |
| { |
| //swap src and dest arrays with particles |
| std::swap(cells, cells2); |
| //swap src and dest arrays with counts of particles |
| std::swap(cnumPars, cnumPars2); |
| // Note, in parallel versions the above swaps may |
| // occure outside RebuildGrid() |
| |
| //initialize destination data structures |
| memset(cnumPars, 0, numCells*sizeof(int)); |
| for(int i=0; i<numCells; i++) |
| { |
| cells[i].next = NULL; |
| last_cells[i] = &cells[i]; |
| } |
| |
| //iterate through source cell lists |
| for(int i = 0; i < numCells; ++i) |
| { |
| Cell *cell2 = &cells2[i]; |
| int np2 = cnumPars2[i]; |
| //iterate through source particles |
| for(int j = 0; j < np2; ++j) |
| { |
| //get destination for source particle |
| int ci = (int)((cell2->p[j % PARTICLES_PER_CELL].x - domainMin.x) / delta.x); |
| int cj = (int)((cell2->p[j % PARTICLES_PER_CELL].y - domainMin.y) / delta.y); |
| int ck = (int)((cell2->p[j % PARTICLES_PER_CELL].z - domainMin.z) / delta.z); |
| // confine to domain |
| // Note, if ProcessCollisions() is working properly these tests are useless |
| if(ci < 0) ci = 0; else if(ci >= nx) ci = nx-1; |
| if(cj < 0) cj = 0; else if(cj >= ny) cj = ny-1; |
| if(ck < 0) ck = 0; else if(ck >= nz) ck = nz-1; |
| |
| #ifdef ENABLE_CFL_CHECK |
| //check that source cell is a neighbor of destination cell |
| bool cfl_cond_satisfied=false; |
| for(int di = -1; di <= 1; ++di) |
| for(int dj = -1; dj <= 1; ++dj) |
| for(int dk = -1; dk <= 1; ++dk) |
| { |
| int ii = ci + di; |
| int jj = cj + dj; |
| int kk = ck + dk; |
| if(ii >= 0 && ii < nx && jj >= 0 && jj < ny && kk >= 0 && kk < nz) |
| { |
| int index = (kk*ny + jj)*nx + ii; |
| if(index == i) |
| { |
| cfl_cond_satisfied=true; |
| break; |
| } |
| } |
| } |
| if(!cfl_cond_satisfied) |
| { |
| std::cerr << "FATAL ERROR: Courant–Friedrichs–Lewy condition not satisfied." << std::endl; |
| exit(1); |
| } |
| #endif //ENABLE_CFL_CHECK |
| |
| //get last pointer in correct destination cell list |
| int index = (ck*ny + cj)*nx + ci; |
| Cell *cell = last_cells[index]; |
| int np = cnumPars[index]; |
| |
| //add another cell structure if everything full |
| if( (np % PARTICLES_PER_CELL == 0) && (cnumPars[index] != 0) ) { |
| cell->next = cellpool_getcell(&pool); |
| cell = cell->next; |
| last_cells[index] = cell; |
| } |
| ++cnumPars[index]; |
| |
| //copy source to destination particle |
| cell->p[np % PARTICLES_PER_CELL].x = cell2->p[j % PARTICLES_PER_CELL].x; |
| cell->p[np % PARTICLES_PER_CELL].y = cell2->p[j % PARTICLES_PER_CELL].y; |
| cell->p[np % PARTICLES_PER_CELL].z = cell2->p[j % PARTICLES_PER_CELL].z; |
| cell->hv[np % PARTICLES_PER_CELL].x = cell2->hv[j % PARTICLES_PER_CELL].x; |
| cell->hv[np % PARTICLES_PER_CELL].y = cell2->hv[j % PARTICLES_PER_CELL].y; |
| cell->hv[np % PARTICLES_PER_CELL].z = cell2->hv[j % PARTICLES_PER_CELL].z; |
| cell->v[np % PARTICLES_PER_CELL].x = cell2->v[j % PARTICLES_PER_CELL].x; |
| cell->v[np % PARTICLES_PER_CELL].y = cell2->v[j % PARTICLES_PER_CELL].y; |
| cell->v[np % PARTICLES_PER_CELL].z = cell2->v[j % PARTICLES_PER_CELL].z; |
| |
| //move pointer to next source cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) { |
| Cell *temp = cell2; |
| cell2 = cell2->next; |
| //return cells to pool that are not statically allocated head of lists |
| if(temp != &cells2[i]) { |
| cellpool_returncell(&pool, temp); |
| } |
| } |
| |
| } // for(int j = 0; j < np2; ++j) |
| //return cells to pool that are not statically allocated head of lists |
| if((cell2 != NULL) && (cell2 != &cells2[i])) { |
| cellpool_returncell(&pool, cell2); |
| } |
| } // for(int i = 0; i < numCells; ++i) |
| } |
| |
| //////////////////////////////////////////////////////////////////////////////// |
| |
| int GetNeighborCells(int ci, int cj, int ck, int *neighCells) |
| { |
| int numNeighCells = 0; |
| |
| // have the nearest particles first -> help branch prediction |
| int my_index = (ck*ny + cj)*nx + ci; |
| neighCells[numNeighCells] = my_index; |
| ++numNeighCells; |
| |
| for(int di = -1; di <= 1; ++di) |
| for(int dj = -1; dj <= 1; ++dj) |
| for(int dk = -1; dk <= 1; ++dk) |
| { |
| int ii = ci + di; |
| int jj = cj + dj; |
| int kk = ck + dk; |
| if(ii >= 0 && ii < nx && jj >= 0 && jj < ny && kk >= 0 && kk < nz) |
| { |
| int index = (kk*ny + jj)*nx + ii; |
| if((index < my_index) && (cnumPars[index] != 0)) |
| { |
| neighCells[numNeighCells] = index; |
| ++numNeighCells; |
| } |
| } |
| } |
| |
| return numNeighCells; |
| } |
| |
| //////////////////////////////////////////////////////////////////////////////// |
| |
| void ComputeForces() |
| { |
| for(int i = 0; i < numCells; ++i) |
| { |
| Cell *cell = &cells[i]; |
| int np = cnumPars[i]; |
| for(int j = 0; j < np; ++j) |
| { |
| cell->density[j % PARTICLES_PER_CELL] = 0.0; |
| cell->a[j % PARTICLES_PER_CELL] = externalAcceleration; |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) { |
| cell = cell->next; |
| } |
| } |
| } |
| |
| int neighCells[3*3*3]; |
| |
| int cindex = 0; |
| for(int ck = 0; ck < nz; ++ck) |
| for(int cj = 0; cj < ny; ++cj) |
| for(int ci = 0; ci < nx; ++ci, ++cindex) |
| { |
| int np = cnumPars[cindex]; |
| if(np == 0) |
| continue; |
| |
| int numNeighCells = GetNeighborCells(ci, cj, ck, neighCells); |
| |
| Cell *cell = &cells[cindex]; |
| for(int ipar = 0; ipar < np; ++ipar) |
| { |
| for(int inc = 0; inc < numNeighCells; ++inc) |
| { |
| int cindexNeigh = neighCells[inc]; |
| Cell *neigh = &cells[cindexNeigh]; |
| int numNeighPars = cnumPars[cindexNeigh]; |
| for(int iparNeigh = 0; iparNeigh < numNeighPars; ++iparNeigh) |
| { |
| //Check address to make sure densities are computed only once per pair |
| if(&neigh->p[iparNeigh % PARTICLES_PER_CELL] < &cell->p[ipar % PARTICLES_PER_CELL]) |
| { |
| fptype distSq = (cell->p[ipar % PARTICLES_PER_CELL] - neigh->p[iparNeigh % PARTICLES_PER_CELL]).GetLengthSq(); |
| if(distSq < hSq) |
| { |
| fptype t = hSq - distSq; |
| fptype tc = t*t*t; |
| cell->density[ipar % PARTICLES_PER_CELL] += tc; |
| neigh->density[iparNeigh % PARTICLES_PER_CELL] += tc; |
| } |
| } |
| //move pointer to next cell in list if end of array is reached |
| if(iparNeigh % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) { |
| neigh = neigh->next; |
| } |
| } |
| } |
| //move pointer to next cell in list if end of array is reached |
| if(ipar % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) { |
| cell = cell->next; |
| } |
| } |
| } |
| |
| const fptype tc = hSq*hSq*hSq; |
| for(int i = 0; i < numCells; ++i) |
| { |
| Cell *cell = &cells[i]; |
| int np = cnumPars[i]; |
| for(int j = 0; j < np; ++j) |
| { |
| cell->density[j % PARTICLES_PER_CELL] += tc; |
| cell->density[j % PARTICLES_PER_CELL] *= densityCoeff; |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) { |
| cell = cell->next; |
| } |
| } |
| } |
| |
| cindex = 0; |
| for(int ck = 0; ck < nz; ++ck) |
| for(int cj = 0; cj < ny; ++cj) |
| for(int ci = 0; ci < nx; ++ci, ++cindex) |
| { |
| int np = cnumPars[cindex]; |
| if(np == 0) |
| continue; |
| |
| int numNeighCells = GetNeighborCells(ci, cj, ck, neighCells); |
| |
| Cell *cell = &cells[cindex]; |
| for(int ipar = 0; ipar < np; ++ipar) |
| { |
| for(int inc = 0; inc < numNeighCells; ++inc) |
| { |
| int cindexNeigh = neighCells[inc]; |
| Cell *neigh = &cells[cindexNeigh]; |
| int numNeighPars = cnumPars[cindexNeigh]; |
| for(int iparNeigh = 0; iparNeigh < numNeighPars; ++iparNeigh) |
| { |
| //Check address to make sure forces are computed only once per pair |
| if(&neigh->p[iparNeigh % PARTICLES_PER_CELL] < &cell->p[ipar % PARTICLES_PER_CELL]) |
| { |
| Vec3 disp = cell->p[ipar % PARTICLES_PER_CELL] - neigh->p[iparNeigh % PARTICLES_PER_CELL]; |
| fptype distSq = disp.GetLengthSq(); |
| if(distSq < hSq) |
| { |
| #ifndef ENABLE_DOUBLE_PRECISION |
| fptype dist = sqrtf(std::max(distSq, (fptype)1e-12)); |
| #else |
| fptype dist = sqrt(std::max(distSq, 1e-12)); |
| #endif //ENABLE_DOUBLE_PRECISION |
| fptype hmr = h - dist; |
| |
| Vec3 acc = disp * pressureCoeff * (hmr*hmr/dist) * (cell->density[ipar % PARTICLES_PER_CELL]+neigh->density[iparNeigh % PARTICLES_PER_CELL] - doubleRestDensity); |
| acc += (neigh->v[iparNeigh % PARTICLES_PER_CELL] - cell->v[ipar % PARTICLES_PER_CELL]) * viscosityCoeff * hmr; |
| acc /= cell->density[ipar % PARTICLES_PER_CELL] * neigh->density[iparNeigh % PARTICLES_PER_CELL]; |
| |
| cell->a[ipar % PARTICLES_PER_CELL] += acc; |
| neigh->a[iparNeigh % PARTICLES_PER_CELL] -= acc; |
| } |
| } |
| //move pointer to next cell in list if end of array is reached |
| if(iparNeigh % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) { |
| neigh = neigh->next; |
| } |
| } |
| } |
| //move pointer to next cell in list if end of array is reached |
| if(ipar % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) { |
| cell = cell->next; |
| } |
| } |
| } |
| } |
| |
| //////////////////////////////////////////////////////////////////////////////// |
| // ProcessCollisions() with container walls |
| // Under the assumptions that |
| // a) a particle will not penetrate a wall |
| // b) a particle will not migrate further than once cell |
| // c) the parSize is smaller than a cell |
| // then only the particles at the perimiters may be influenced by the walls |
| #if 0 |
| void ProcessCollisions() |
| { |
| for(int i = 0; i < numCells; ++i) |
| { |
| Cell *cell = &cells[i]; |
| int np = cnumPars[i]; |
| for(int j = 0; j < np; ++j) |
| { |
| Vec3 pos = cell->p[j % PARTICLES_PER_CELL] + cell->hv[j % PARTICLES_PER_CELL] * timeStep; |
| |
| fptype diff = parSize - (pos.x - domainMin.x); |
| if(diff > epsilon) |
| cell->a[j % PARTICLES_PER_CELL].x += stiffnessCollisions*diff - damping*cell->v[j % PARTICLES_PER_CELL].x; |
| |
| diff = parSize - (domainMax.x - pos.x); |
| if(diff > epsilon) |
| cell->a[j % PARTICLES_PER_CELL].x -= stiffnessCollisions*diff + damping*cell->v[j % PARTICLES_PER_CELL].x; |
| |
| diff = parSize - (pos.y - domainMin.y); |
| if(diff > epsilon) |
| cell->a[j % PARTICLES_PER_CELL].y += stiffnessCollisions*diff - damping*cell->v[j % PARTICLES_PER_CELL].y; |
| |
| diff = parSize - (domainMax.y - pos.y); |
| if(diff > epsilon) |
| cell->a[j % PARTICLES_PER_CELL].y -= stiffnessCollisions*diff + damping*cell->v[j % PARTICLES_PER_CELL].y; |
| |
| diff = parSize - (pos.z - domainMin.z); |
| if(diff > epsilon) |
| cell->a[j % PARTICLES_PER_CELL].z += stiffnessCollisions*diff - damping*cell->v[j % PARTICLES_PER_CELL].z; |
| |
| diff = parSize - (domainMax.z - pos.z); |
| if(diff > epsilon) |
| cell->a[j % PARTICLES_PER_CELL].z -= stiffnessCollisions*diff + damping*cell->v[j % PARTICLES_PER_CELL].z; |
| |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) |
| cell = cell->next; |
| } |
| } |
| } |
| #else |
| // Notes on USE_ImpeneratableWall |
| // When particle is detected beyond cell wall it is repositioned at cell wall |
| // velocity is not changed, thus conserving momentum. |
| // What this means though it the prior AdvanceParticles had positioned the |
| // particle beyond the cell wall and thus the visualization will show these |
| // as artifacts. The proper place for USE_ImpeneratableWall is after AdvanceParticles. |
| // This would entail a 2nd pass on the perimiters after AdvanceParticles (as opposed |
| // to inside AdvanceParticles). Your fluid dynamisist should properly devise the |
| // equasions. |
| void ProcessCollisions() |
| { |
| int x,y,z, index; |
| x=0; // along the domainMin.x wall |
| for(y=0; y<ny; ++y) |
| { |
| for(z=0; z<nz; ++z) |
| { |
| int ci = (z*ny + y)*nx + x; |
| Cell *cell = &cells[ci]; |
| int np = cnumPars[ci]; |
| for(int j = 0; j < np; ++j) |
| { |
| int ji = j % PARTICLES_PER_CELL; |
| fptype pos_x = cell->p[ji].x + cell->hv[ji].x * timeStep; |
| fptype diff = parSize - (pos_x - domainMin.x); |
| if(diff > epsilon) |
| cell->a[ji].x += stiffnessCollisions*diff - damping*cell->v[ji].x; |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) |
| cell = cell->next; |
| } |
| } |
| } |
| x=nx-1; // along the domainMax.x wall |
| for(y=0; y<ny; ++y) |
| { |
| for(z=0; z<nz; ++z) |
| { |
| int ci = (z*ny + y)*nx + x; |
| Cell *cell = &cells[ci]; |
| int np = cnumPars[ci]; |
| for(int j = 0; j < np; ++j) |
| { |
| int ji = j % PARTICLES_PER_CELL; |
| fptype pos_x = cell->p[ji].x + cell->hv[ji].x * timeStep; |
| fptype diff = parSize - (domainMax.x - pos_x); |
| if(diff > epsilon) |
| cell->a[ji].x -= stiffnessCollisions*diff + damping*cell->v[ji].x; |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) |
| cell = cell->next; |
| } |
| } |
| } |
| y=0; // along the domainMin.y wall |
| for(x=0; x<nx; ++x) |
| { |
| for(z=0; z<nz; ++z) |
| { |
| int ci = (z*ny + y)*nx + x; |
| Cell *cell = &cells[ci]; |
| int np = cnumPars[ci]; |
| for(int j = 0; j < np; ++j) |
| { |
| int ji = j % PARTICLES_PER_CELL; |
| fptype pos_y = cell->p[ji].y + cell->hv[ji].y * timeStep; |
| fptype diff = parSize - (pos_y - domainMin.y); |
| if(diff > epsilon) |
| cell->a[ji].y += stiffnessCollisions*diff - damping*cell->v[ji].y; |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) |
| cell = cell->next; |
| } |
| } |
| } |
| y=ny-1; // along the domainMax.y wall |
| for(x=0; x<nx; ++x) |
| { |
| for(z=0; z<nz; ++z) |
| { |
| int ci = (z*ny + y)*nx + x; |
| Cell *cell = &cells[ci]; |
| int np = cnumPars[ci]; |
| for(int j = 0; j < np; ++j) |
| { |
| int ji = j % PARTICLES_PER_CELL; |
| fptype pos_y = cell->p[ji].y + cell->hv[ji].y * timeStep; |
| fptype diff = parSize - (domainMax.y - pos_y); |
| if(diff > epsilon) |
| cell->a[ji].y -= stiffnessCollisions*diff + damping*cell->v[ji].y; |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) |
| cell = cell->next; |
| } |
| } |
| } |
| z=0; // along the domainMin.z wall |
| for(x=0; x<nx; ++x) |
| { |
| for(y=0; y<ny; ++y) |
| { |
| int ci = (z*ny + y)*nx + x; |
| Cell *cell = &cells[ci]; |
| int np = cnumPars[ci]; |
| for(int j = 0; j < np; ++j) |
| { |
| int ji = j % PARTICLES_PER_CELL; |
| fptype pos_z = cell->p[ji].z + cell->hv[ji].z * timeStep; |
| fptype diff = parSize - (pos_z - domainMin.z); |
| if(diff > epsilon) |
| cell->a[ji].z += stiffnessCollisions*diff - damping*cell->v[ji].z; |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) |
| cell = cell->next; |
| } |
| } |
| } |
| z=nz-1; // along the domainMax.z wall |
| for(x=0; x<nx; ++x) |
| { |
| for(y=0; y<ny; ++y) |
| { |
| int ci = (z*ny + y)*nx + x; |
| Cell *cell = &cells[ci]; |
| int np = cnumPars[ci]; |
| for(int j = 0; j < np; ++j) |
| { |
| int ji = j % PARTICLES_PER_CELL; |
| fptype pos_z = cell->p[ji].z + cell->hv[ji].z * timeStep; |
| fptype diff = parSize - (domainMax.z - pos_z); |
| if(diff > epsilon) |
| cell->a[ji].z -= stiffnessCollisions*diff + damping*cell->v[ji].z; |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) |
| cell = cell->next; |
| } |
| } |
| } |
| } |
| #define USE_ImpeneratableWall |
| #if defined(USE_ImpeneratableWall) |
| void ProcessCollisions2() |
| { |
| int x,y,z, index; |
| x=0; // along the domainMin.x wall |
| for(y=0; y<ny; ++y) |
| { |
| for(z=0; z<nz; ++z) |
| { |
| int ci = (z*ny + y)*nx + x; |
| Cell *cell = &cells[ci]; |
| int np = cnumPars[ci]; |
| for(int j = 0; j < np; ++j) |
| { |
| int ji = j % PARTICLES_PER_CELL; |
| fptype diff = cell->p[ji].x - domainMin.x; |
| if(diff < Zero) |
| { |
| cell->p[ji].x = domainMin.x - diff; |
| cell->v[ji].x = -cell->v[ji].x; |
| cell->hv[ji].x = -cell->hv[ji].x; |
| } |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) |
| cell = cell->next; |
| } |
| } |
| } |
| x=nx-1; // along the domainMax.x wall |
| for(y=0; y<ny; ++y) |
| { |
| for(z=0; z<nz; ++z) |
| { |
| int ci = (z*ny + y)*nx + x; |
| Cell *cell = &cells[ci]; |
| int np = cnumPars[ci]; |
| for(int j = 0; j < np; ++j) |
| { |
| int ji = j % PARTICLES_PER_CELL; |
| fptype diff = domainMax.x - cell->p[ji].x; |
| if(diff < Zero) |
| { |
| cell->p[ji].x = domainMax.x + diff; |
| cell->v[ji].x = -cell->v[ji].x; |
| cell->hv[ji].x = -cell->hv[ji].x; |
| } |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) |
| cell = cell->next; |
| } |
| } |
| } |
| y=0; // along the domainMin.y wall |
| for(x=0; x<nx; ++x) |
| { |
| for(z=0; z<nz; ++z) |
| { |
| int ci = (z*ny + y)*nx + x; |
| Cell *cell = &cells[ci]; |
| int np = cnumPars[ci]; |
| for(int j = 0; j < np; ++j) |
| { |
| int ji = j % PARTICLES_PER_CELL; |
| fptype diff = cell->p[ji].y - domainMin.y; |
| if(diff < Zero) |
| { |
| cell->p[ji].y = domainMin.y - diff; |
| cell->v[ji].y = -cell->v[ji].y; |
| cell->hv[ji].y = -cell->hv[ji].y; |
| } |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) |
| cell = cell->next; |
| } |
| } |
| } |
| y=ny-1; // along the domainMax.y wall |
| for(x=0; x<nx; ++x) |
| { |
| for(z=0; z<nz; ++z) |
| { |
| int ci = (z*ny + y)*nx + x; |
| Cell *cell = &cells[ci]; |
| int np = cnumPars[ci]; |
| for(int j = 0; j < np; ++j) |
| { |
| int ji = j % PARTICLES_PER_CELL; |
| fptype diff = domainMax.y - cell->p[ji].y; |
| if(diff < Zero) |
| { |
| cell->p[ji].y = domainMax.y + diff; |
| cell->v[ji].y = -cell->v[ji].y; |
| cell->hv[ji].y = -cell->hv[ji].y; |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) |
| cell = cell->next; |
| } |
| } |
| } |
| } |
| z=0; // along the domainMin.z wall |
| for(x=0; x<nx; ++x) |
| { |
| for(y=0; y<ny; ++y) |
| { |
| int ci = (z*ny + y)*nx + x; |
| Cell *cell = &cells[ci]; |
| int np = cnumPars[ci]; |
| for(int j = 0; j < np; ++j) |
| { |
| int ji = j % PARTICLES_PER_CELL; |
| fptype diff = cell->p[ji].z - domainMin.z; |
| if(diff < Zero) |
| { |
| cell->p[ji].z = domainMin.z - diff; |
| cell->v[ji].z = -cell->v[ji].z; |
| cell->hv[ji].z = -cell->hv[ji].z; |
| } |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) |
| cell = cell->next; |
| } |
| } |
| } |
| z=nz-1; // along the domainMax.z wall |
| for(x=0; x<nx; ++x) |
| { |
| for(y=0; y<ny; ++y) |
| { |
| int ci = (z*ny + y)*nx + x; |
| Cell *cell = &cells[ci]; |
| int np = cnumPars[ci]; |
| for(int j = 0; j < np; ++j) |
| { |
| int ji = j % PARTICLES_PER_CELL; |
| fptype diff = domainMax.z - cell->p[ji].z; |
| if(diff < Zero) |
| { |
| cell->p[ji].z = domainMax.z + diff; |
| cell->v[ji].z = -cell->v[ji].z; |
| cell->hv[ji].z = -cell->hv[ji].z; |
| } |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) |
| cell = cell->next; |
| } |
| } |
| } |
| } |
| #endif |
| #endif |
| //////////////////////////////////////////////////////////////////////////////// |
| |
| void AdvanceParticles() |
| { |
| for(int i = 0; i < numCells; ++i) |
| { |
| Cell *cell = &cells[i]; |
| int np = cnumPars[i]; |
| for(int j = 0; j < np; ++j) |
| { |
| Vec3 v_half = cell->hv[j % PARTICLES_PER_CELL] + cell->a[j % PARTICLES_PER_CELL]*timeStep; |
| #if defined(USE_ImpeneratableWall) |
| // N.B. The integration of the position can place the particle |
| // outside the domain. Although we could place a test in this loop |
| // we would be unnecessarily testing particles on interior cells. |
| // Therefore, to reduce the amount of computations we make a later |
| // pass on the perimiter cells to account for particle migration |
| // beyond domain |
| #endif |
| cell->p[j % PARTICLES_PER_CELL] += v_half * timeStep; |
| cell->v[j % PARTICLES_PER_CELL] = cell->hv[j % PARTICLES_PER_CELL] + v_half; |
| cell->v[j % PARTICLES_PER_CELL] *= 0.5; |
| cell->hv[j % PARTICLES_PER_CELL] = v_half; |
| |
| //move pointer to next cell in list if end of array is reached |
| if(j % PARTICLES_PER_CELL == PARTICLES_PER_CELL-1) { |
| cell = cell->next; |
| } |
| } |
| } |
| } |
| |
| //////////////////////////////////////////////////////////////////////////////// |
| |
| void AdvanceFrame() |
| { |
| RebuildGrid(); |
| ComputeForces(); |
| ProcessCollisions(); |
| AdvanceParticles(); |
| #if defined(USE_ImpeneratableWall) |
| // N.B. The integration of the position can place the particle |
| // outside the domain. We now make a pass on the perimiter cells |
| // to account for particle migration beyond domain. |
| ProcessCollisions2(); |
| #endif |
| #ifdef ENABLE_STATISTICS |
| float mean, stddev; |
| int i; |
| |
| mean = (float)numParticles/(float)numCells; |
| stddev = 0.0; |
| for(i=0; i<numCells; i++) { |
| stddev += (mean-cnumPars[i])*(mean-cnumPars[i]); |
| } |
| stddev = sqrtf(stddev); |
| std::cout << "Cell statistics: mean=" << mean << " particles, stddev=" << stddev << " particles." << std::endl; |
| #endif |
| } |
| |
| //////////////////////////////////////////////////////////////////////////////// |
| |
| int main(int argc, char *argv[]) |
| { |
| #ifdef PARSEC_VERSION |
| #define __PARSEC_STRING(x) #x |
| #define __PARSEC_XSTRING(x) __PARSEC_STRING(x) |
| std::cout << "PARSEC Benchmark Suite Version "__PARSEC_XSTRING(PARSEC_VERSION) << std::endl << std::flush; |
| #else |
| std::cout << "PARSEC Benchmark Suite" << std::endl << std::flush; |
| #endif //PARSEC_VERSION |
| #ifdef ENABLE_PARSEC_HOOKS |
| __parsec_bench_begin(__parsec_fluidanimate); |
| #endif |
| |
| if(argc < 4 || argc >= 6) |
| { |
| std::cout << "Usage: " << argv[0] << " <threadnum> <framenum> <.fluid input file> [.fluid output file]" << std::endl; |
| return -1; |
| } |
| |
| int threadnum = atoi(argv[1]); |
| int framenum = atoi(argv[2]); |
| |
| //Check arguments |
| if(threadnum != 1) { |
| std::cerr << "<threadnum> must be 1 (serial version)" << std::endl; |
| return -1; |
| } |
| if(framenum < 1) { |
| std::cerr << "<framenum> must at least be 1" << std::endl; |
| return -1; |
| } |
| |
| #ifdef ENABLE_CFL_CHECK |
| std::cout << "WARNING: Check for Courant–Friedrichs–Lewy condition enabled. Do not use for performance measurements." << std::endl; |
| #endif |
| |
| InitSim(argv[3]); |
| #ifdef ENABLE_VISUALIZATION |
| InitVisualizationMode(&argc, argv, &AdvanceFrame, &numCells, &cells, &cnumPars); |
| #endif |
| |
| #ifndef ENABLE_VISUALIZATION |
| |
| //core of benchmark program (the Region-of-Interest) |
| #ifdef ENABLE_PARSEC_HOOKS |
| __parsec_roi_begin(); |
| #endif |
| for(int i = 0; i < framenum; ++i) |
| AdvanceFrame(); |
| #ifdef ENABLE_PARSEC_HOOKS |
| __parsec_roi_end(); |
| #endif |
| |
| #else //ENABLE_VISUALIZATION |
| Visualize(); |
| #endif //ENABLE_VISUALIZATION |
| |
| if(argc > 4) |
| SaveFile(argv[4]); |
| CleanUpSim(); |
| |
| #ifdef ENABLE_PARSEC_HOOKS |
| __parsec_bench_end(); |
| #endif |
| |
| return 0; |
| } |
| |
| //////////////////////////////////////////////////////////////////////////////// |
| |
