intern/smoke/intern/WTURBULENCE.cpp

   1 /** \file smoke/intern/WTURBULENCE.cpp
   2  *  \ingroup smoke
   3  */
   4 //////////////////////////////////////////////////////////////////////
   5 // This file is part of Wavelet Turbulence.
   6 //
   7 // Wavelet Turbulence is free software: you can redistribute it and/or modify
   8 // it under the terms of the GNU General Public License as published by
   9 // the Free Software Foundation, either version 3 of the License, or
  10 // (at your option) any later version.
  11 //
  12 // Wavelet Turbulence is distributed in the hope that it will be useful,
  13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 // GNU General Public License for more details.
  16 //
  17 // You should have received a copy of the GNU General Public License
  18 // along with Wavelet Turbulence.  If not, see <http://www.gnu.org/licenses/>.
  19 //
  20 // Copyright 2008 Theodore Kim and Nils Thuerey
  21 //
  22 // WTURBULENCE handling
  23 ///////////////////////////////////////////////////////////////////////////////////
  24 // Parallelized turbulence even further. TNT matrix library functions
  25 // rewritten to improve performance.
  26 //              - MiikaH
  27 //////////////////////////////////////////////////////////////////////
  28
  29 #include "WTURBULENCE.h"
  30 #include "INTERPOLATE.h"
  31 #include "IMAGE.h"
  32 #include <MERSENNETWISTER.h>
  33 #include "WAVELET_NOISE.h"
  34 #include "FFT_NOISE.h"
  35 #include "EIGENVALUE_HELPER.h"
  36 #include "LU_HELPER.h"
  37 #include "SPHERE.h"
  38 #include <zlib.h>
  39 #include <math.h>
  40
  41 // needed to access static advection functions
  42 #include "FLUID_3D.h"
  43
  44 #if PARALLEL==1
  45 #include <omp.h>
  46 #endif // PARALLEL
  47
  48 // 2^ {-5/6}
  49 static const float persistence = 0.56123f;
  50
  51 //////////////////////////////////////////////////////////////////////
  52 // constructor
  53 //////////////////////////////////////////////////////////////////////
  54 WTURBULENCE::WTURBULENCE(int xResSm, int yResSm, int zResSm, int amplify, int noisetype)
  55 {
  56         // if noise magnitude is below this threshold, its contribution
  57         // is negilgible, so stop evaluating new octaves
  58         _cullingThreshold = 1e-3;
  59
  60         // factor by which to increase the simulation resolution
  61         _amplify = amplify;
  62
  63         // manually adjust the overall amount of turbulence
  64         // DG - RNA-fied _strength = 2.;
  65
  66         // add the corresponding octaves of noise
  67         _octaves = (int)(log((float)_amplify) / log(2.0f) + 0.5); // XXX DEBUG/ TODO: int casting correct? - dg
  68
  69         // noise resolution
  70         _xResBig = _amplify * xResSm;
  71         _yResBig = _amplify * yResSm;
  72         _zResBig = _amplify * zResSm;
  73         _resBig = Vec3Int(_xResBig, _yResBig, _zResBig);
  74         _invResBig = Vec3(1./(float)_resBig[0], 1./(float)_resBig[1], 1./(float)_resBig[2]);
  75         _slabSizeBig = _xResBig*_yResBig;
  76         _totalCellsBig = _slabSizeBig * _zResBig;
  77
  78         // original / small resolution
  79         _xResSm = xResSm;
  80         _yResSm = yResSm;
  81         _zResSm = zResSm;
  82         _resSm = Vec3Int(xResSm, yResSm, zResSm);
  83         _invResSm = Vec3(1./(float)_resSm[0], 1./(float)_resSm[1], 1./(float)_resSm[2] );
  84         _slabSizeSm = _xResSm*_yResSm;
  85         _totalCellsSm = _slabSizeSm * _zResSm;
  86
  87         // allocate high resolution density field
  88         _totalStepsBig = 0;
  89         _densityBig = new float[_totalCellsBig];
  90         _densityBigOld = new float[_totalCellsBig];
  91
  92         for(int i = 0; i < _totalCellsBig; i++) {
  93                 _densityBig[i] =
  94                 _densityBigOld[i] = 0.;
  95         }
  96
  97         // allocate & init texture coordinates
  98         _tcU = new float[_totalCellsSm];
  99         _tcV = new float[_totalCellsSm];
 100         _tcW = new float[_totalCellsSm];
 101         _tcTemp = new float[_totalCellsSm];
 102
 103         // map all
 104         const float dx = 1./(float)(_resSm[0]);
 105         const float dy = 1./(float)(_resSm[1]);
 106         const float dz = 1./(float)(_resSm[2]);
 107         int index = 0;
 108         for (int z = 0; z < _zResSm; z++)
 109         for (int y = 0; y < _yResSm; y++)
 110                 for (int x = 0; x < _xResSm; x++, index++)
 111                 {
 112                 _tcU[index] = x*dx;
 113                 _tcV[index] = y*dy;
 114                 _tcW[index] = z*dz;
 115                 _tcTemp[index] = 0.;
 116                 }
 117
 118         // noise tiles
 119         _noiseTile = new float[noiseTileSize * noiseTileSize * noiseTileSize];
 120         /*
 121         std::string noiseTileFilename = std::string("noise.wavelets");
 122         generateTile_WAVELET(_noiseTile, noiseTileFilename);
 123         */
 124         setNoise(noisetype);
 125         /*
 126         std::string noiseTileFilename = std::string("noise.fft");
 127         generatTile_FFT(_noiseTile, noiseTileFilename);
 128         */
 129 }
 130
 131 //////////////////////////////////////////////////////////////////////
 132 // destructor
 133 //////////////////////////////////////////////////////////////////////
 134 WTURBULENCE::~WTURBULENCE() {
 135   delete[] _densityBig;
 136   delete[] _densityBigOld;
 137
 138   delete[] _tcU;
 139   delete[] _tcV;
 140   delete[] _tcW;
 141   delete[] _tcTemp;
 142
 143   delete[] _noiseTile;
 144 }
 145
 146 //////////////////////////////////////////////////////////////////////
 147 // Change noise type
 148 //
 149 // type (1<<0) = wavelet / 2
 150 // type (1<<1) = FFT / 4
 151 // type (1<<2) = curl / 8
 152 //////////////////////////////////////////////////////////////////////
 153 void WTURBULENCE::setNoise(int type)
 154 {
 155         if(type == (1<<1)) // FFT
 156         {
 157                 // needs fft
 158                 #if FFTW3==1
 159                 std::string noiseTileFilename = std::string("noise.fft");
 160                 generatTile_FFT(_noiseTile, noiseTileFilename);
 161                 #endif
 162         }
 163         else if(type == (1<<2)) // curl
 164         {
 165                 // TODO: not supported yet
 166         }
 167         else // standard - wavelet
 168         {
 169                 std::string noiseTileFilename = std::string("noise.wavelets");
 170                 generateTile_WAVELET(_noiseTile, noiseTileFilename);
 171         }
 172 }
 173
 174 // init direct access functions from blender
 175 void WTURBULENCE::initBlenderRNA(float *strength)
 176 {
 177         _strength = strength;
 178 }
 179
 180 //////////////////////////////////////////////////////////////////////
 181 // Get the smallest valid x derivative
 182 //
 183 // Takes the one-sided finite difference in both directions and
 184 // selects the smaller of the two
 185 //////////////////////////////////////////////////////////////////////
 186 static float minDx(int x, int y, int z, float* input, Vec3Int res)
 187 {
 188   const int index = x + y * res[0] + z * res[0] * res[1];
 189   const int maxx = res[0]-2;
 190
 191   // get grid values
 192   float center = input[index];
 193   float left  = (x <= 1)    ? FLT_MAX : input[index - 1];
 194   float right = (x >= maxx) ? FLT_MAX : input[index + 1];
 195
 196   const float dx = res[0];
 197
 198   // get all the derivative estimates
 199   float dLeft   = (x <= 1)     ? FLT_MAX : (center - left) * dx;
 200   float dRight  = (x >= maxx)  ? FLT_MAX : (right - center) * dx;
 201   float dCenter = (x <= 1 || x >= maxx) ? FLT_MAX : (right - left) * dx * 0.5f;
 202
 203   // if it's on a boundary, only one estimate is valid
 204   if (x <= 1) return dRight;
 205   if (x >= maxx) return dLeft;
 206
 207   // if it's not on a boundary, get the smallest one
 208   float finalD;
 209   finalD = (fabs(dCenter) < fabs(dRight)) ? dCenter : dRight;
 210   finalD = (fabs(finalD)  < fabs(dLeft))  ? finalD  : dLeft;
 211
 212   return finalD;
 213 }
 214
 215 //////////////////////////////////////////////////////////////////////
 216 // get the smallest valid y derivative
 217 //
 218 // Takes the one-sided finite difference in both directions and
 219 // selects the smaller of the two
 220 //////////////////////////////////////////////////////////////////////
 221 static float minDy(int x, int y, int z, float* input, Vec3Int res)
 222 {
 223   const int index = x + y * res[0] + z * res[0] * res[1];
 224   const int maxy = res[1]-2;
 225
 226   // get grid values
 227   float center = input[index];
 228   float down  = (y <= 1) ? FLT_MAX : input[index - res[0]];
 229   float up = (y >= maxy) ? FLT_MAX : input[index + res[0]];
 230
 231   const float dx = res[1]; // only for square domains
 232
 233   // get all the derivative estimates
 234   float dDown   = (y <= 1)  ? FLT_MAX : (center - down) * dx;
 235   float dUp  = (y >= maxy)  ? FLT_MAX : (up - center) * dx;
 236   float dCenter = (y <= 1 || y >= maxy) ? FLT_MAX : (up - down) * dx * 0.5f;
 237
 238   // if it's on a boundary, only one estimate is valid
 239   if (y <= 1) return dUp;
 240   if (y >= maxy) return dDown;
 241
 242   // if it's not on a boundary, get the smallest one
 243   float finalD = (fabs(dCenter) < fabs(dUp)) ? dCenter : dUp;
 244   finalD = (fabs(finalD) < fabs(dDown)) ? finalD : dDown;
 245
 246   return finalD;
 247 }
 248
 249 //////////////////////////////////////////////////////////////////////
 250 // get the smallest valid z derivative
 251 //
 252 // Takes the one-sided finite difference in both directions and
 253 // selects the smaller of the two
 254 //////////////////////////////////////////////////////////////////////
 255 static float minDz(int x, int y, int z, float* input, Vec3Int res)
 256 {
 257   const int slab = res[0]*res[1];
 258   const int index = x + y * res[0] + z * slab;
 259   const int maxz = res[2]-2;
 260
 261   // get grid values
 262   float center = input[index];
 263   float front  = (z <= 1) ? FLT_MAX : input[index - slab];
 264   float back = (z >= maxz) ? FLT_MAX : input[index + slab];
 265
 266   const float dx = res[2]; // only for square domains
 267
 268   // get all the derivative estimates
 269   float dfront   = (z <= 1)  ? FLT_MAX : (center - front) * dx;
 270   float dback  = (z >= maxz)  ? FLT_MAX : (back - center) * dx;
 271   float dCenter = (z <= 1 || z >= maxz) ? FLT_MAX : (back - front) * dx * 0.5f;
 272
 273   // if it's on a boundary, only one estimate is valid
 274   if (z <= 1) return dback;
 275   if (z >= maxz) return dfront;
 276
 277   // if it's not on a boundary, get the smallest one
 278   float finalD = (fabs(dCenter) < fabs(dback)) ? dCenter : dback;
 279   finalD = (fabs(finalD) < fabs(dfront)) ? finalD : dfront;
 280
 281   return finalD;
 282 }
 283
 284 //////////////////////////////////////////////////////////////////////
 285 // handle texture coordinates (advection, reset, eigenvalues),
 286 // Beware -- uses big density maccormack as temporary arrays
 287 //////////////////////////////////////////////////////////////////////
 288 void WTURBULENCE::advectTextureCoordinates (float dtOrg, float* xvel, float* yvel, float* zvel, float *tempBig1, float *tempBig2) {
 289
 290   // advection
 291   SWAP_POINTERS(_tcTemp, _tcU);
 292   FLUID_3D::copyBorderX(_tcTemp, _resSm, 0 , _resSm[2]);
 293   FLUID_3D::copyBorderY(_tcTemp, _resSm, 0 , _resSm[2]);
 294   FLUID_3D::copyBorderZ(_tcTemp, _resSm, 0 , _resSm[2]);
 295   FLUID_3D::advectFieldMacCormack1(dtOrg, xvel, yvel, zvel,
 296       _tcTemp, tempBig1, _resSm, 0 , _resSm[2]);
 297   FLUID_3D::advectFieldMacCormack2(dtOrg, xvel, yvel, zvel,
 298       _tcTemp, _tcU, tempBig1, tempBig2, _resSm, NULL, 0 , _resSm[2]);
 299
 300   SWAP_POINTERS(_tcTemp, _tcV);
 301   FLUID_3D::copyBorderX(_tcTemp, _resSm, 0 , _resSm[2]);
 302   FLUID_3D::copyBorderY(_tcTemp, _resSm, 0 , _resSm[2]);
 303   FLUID_3D::copyBorderZ(_tcTemp, _resSm, 0 , _resSm[2]);
 304   FLUID_3D::advectFieldMacCormack1(dtOrg, xvel, yvel, zvel,
 305       _tcTemp, tempBig1, _resSm, 0 , _resSm[2]);
 306   FLUID_3D::advectFieldMacCormack2(dtOrg, xvel, yvel, zvel,
 307       _tcTemp, _tcV, tempBig1, tempBig2, _resSm, NULL, 0 , _resSm[2]);
 308
 309   SWAP_POINTERS(_tcTemp, _tcW);
 310   FLUID_3D::copyBorderX(_tcTemp, _resSm, 0 , _resSm[2]);
 311   FLUID_3D::copyBorderY(_tcTemp, _resSm, 0 , _resSm[2]);
 312   FLUID_3D::copyBorderZ(_tcTemp, _resSm, 0 , _resSm[2]);
 313   FLUID_3D::advectFieldMacCormack1(dtOrg, xvel, yvel, zvel,
 314       _tcTemp, tempBig1, _resSm, 0 , _resSm[2]);
 315   FLUID_3D::advectFieldMacCormack2(dtOrg, xvel, yvel, zvel,
 316       _tcTemp, _tcW, tempBig1, tempBig2, _resSm, NULL, 0 , _resSm[2]);
 317 }
 318
 319 //////////////////////////////////////////////////////////////////////
 320 // Compute the eigenvalues of the advected texture
 321 //////////////////////////////////////////////////////////////////////
 322 void WTURBULENCE::computeEigenvalues(float *_eigMin, float *_eigMax) {
 323   // stats
 324   float maxeig = -1.;
 325   float mineig = 10.;
 326
 327   // texture coordinate eigenvalues
 328   for (int z = 1; z < _zResSm-1; z++) {
 329     for (int y = 1; y < _yResSm-1; y++)
 330       for (int x = 1; x < _xResSm-1; x++)
 331       {
 332         const int index = x+ y *_resSm[0] + z*_slabSizeSm;
 333
 334         // compute jacobian
 335         float jacobian[3][3] = {
 336           { minDx(x, y, z, _tcU, _resSm), minDx(x, y, z, _tcV, _resSm), minDx(x, y, z, _tcW, _resSm) } ,
 337           { minDy(x, y, z, _tcU, _resSm), minDy(x, y, z, _tcV, _resSm), minDy(x, y, z, _tcW, _resSm) } ,
 338           { minDz(x, y, z, _tcU, _resSm), minDz(x, y, z, _tcV, _resSm), minDz(x, y, z, _tcW, _resSm) }
 339         };
 340
 341         // ONLY compute the eigenvalues after checking that the matrix
 342         // is nonsingular
 343         sLU LU = computeLU(jacobian);
 344
 345         if (isNonsingular(LU))
 346         {
 347           // get the analytic eigenvalues, quite slow right now...
 348           Vec3 eigenvalues = Vec3(1.);
 349           computeEigenvalues3x3( &eigenvalues[0], jacobian);
 350           _eigMax[index] = MAX3V(eigenvalues);
 351           _eigMin[index] = MIN3V(eigenvalues);
 352           maxeig = MAX(_eigMax[index],maxeig);
 353           mineig = MIN(_eigMin[index],mineig);
 354         }
 355         else
 356         {
 357           _eigMax[index] = 10.0f;
 358           _eigMin[index] = 0.1;
 359         }
 360       }
 361   }
 362 }
 363
 364 //////////////////////////////////////////////////////////////////////
 365 // advect & reset texture coordinates based on eigenvalues
 366 //////////////////////////////////////////////////////////////////////
 367 void WTURBULENCE::resetTextureCoordinates(float *_eigMin, float *_eigMax)
 368 {
 369   // allowed deformation of the textures
 370   const float limit = 2.f;
 371   const float limitInv = 1./limit;
 372
 373   // standard reset
 374   int resets = 0;
 375   const float dx = 1./(float)(_resSm[0]);
 376   const float dy = 1./(float)(_resSm[1]);
 377   const float dz = 1./(float)(_resSm[2]);
 378
 379   for (int z = 1; z < _zResSm-1; z++)
 380     for (int y = 1; y < _yResSm-1; y++)
 381       for (int x = 1; x < _xResSm-1; x++)
 382       {
 383         const int index = x+ y *_resSm[0] + z*_slabSizeSm;
 384         if (_eigMax[index] > limit || _eigMin[index] < limitInv)
 385         {
 386           _tcU[index] = (float)x * dx;
 387           _tcV[index] = (float)y * dy;
 388           _tcW[index] = (float)z * dz;
 389           resets++;
 390         }
 391       }
 392 }
 393
 394 //////////////////////////////////////////////////////////////////////
 395 // Compute the highest frequency component of the wavelet
 396 // decomposition
 397 //////////////////////////////////////////////////////////////////////
 398 void WTURBULENCE::decomposeEnergy(float *_energy, float *_highFreqEnergy)
 399 {
 400   // do the decomposition -- the goal here is to have
 401   // the energy with the high frequency component stomped out
 402   // stored in _tcTemp when it is done. _highFreqEnergy is only used
 403   // as an additional temp array
 404
 405   // downsample input
 406   downsampleXNeumann(_highFreqEnergy, _energy, _xResSm, _yResSm, _zResSm);
 407   downsampleYNeumann(_tcTemp, _highFreqEnergy, _xResSm, _yResSm, _zResSm);
 408   downsampleZNeumann(_highFreqEnergy, _tcTemp, _xResSm, _yResSm, _zResSm);
 409
 410   // upsample input
 411   upsampleZNeumann(_tcTemp, _highFreqEnergy, _xResSm, _yResSm, _zResSm);
 412   upsampleYNeumann(_highFreqEnergy, _tcTemp, _xResSm, _yResSm, _zResSm);
 413   upsampleXNeumann(_tcTemp, _highFreqEnergy, _xResSm, _yResSm, _zResSm);
 414
 415   // subtract the down and upsampled field from the original field --
 416   // what should be left over is solely the high frequency component
 417         int index = 0;
 418   for (int z = 0; z < _zResSm; z++)
 419     for (int y = 0; y < _yResSm; y++) {
 420       for (int x = 0; x < _xResSm; x++, index++) {
 421         // brute force reset of boundaries
 422         if(z >= _zResSm - 1 || x >= _xResSm - 1 || y >= _yResSm - 1 || z <= 0 || y <= 0 || x <= 0)
 423           _highFreqEnergy[index] = 0.;
 424         else
 425           _highFreqEnergy[index] = _energy[index] - _tcTemp[index];
 426     }
 427   }
 428 }
 429
 430 //////////////////////////////////////////////////////////////////////
 431 // compute velocity from energies and march into obstacles
 432 // for wavelet decomposition
 433 //////////////////////////////////////////////////////////////////////
 434 void WTURBULENCE::computeEnergy(float *_energy, float* xvel, float* yvel, float* zvel, unsigned char *obstacles)
 435 {
 436   // compute everywhere
 437   for (int x = 0; x < _totalCellsSm; x++)
 438     _energy[x] = 0.5f * (xvel[x] * xvel[x] + yvel[x] * yvel[x] + zvel[x] * zvel[x]);
 439
 440   FLUID_3D::copyBorderX(_energy, _resSm, 0 , _resSm[2]);
 441   FLUID_3D::copyBorderY(_energy, _resSm, 0 , _resSm[2]);
 442   FLUID_3D::copyBorderZ(_energy, _resSm, 0 , _resSm[2]);
 443
 444   // pseudo-march the values into the obstacles
 445   // the wavelet upsampler only uses a 3x3 support neighborhood, so
 446   // propagating the values in by 4 should be sufficient
 447   int index;
 448
 449   // iterate
 450   for (int iter = 0; iter < 4; iter++)
 451   {
 452     index = _slabSizeSm + _xResSm + 1;
 453     for (int z = 1; z < _zResSm - 1; z++, index += 2 * _xResSm)
 454       for (int y = 1; y < _yResSm - 1; y++, index += 2)
 455         for (int x = 1; x < _xResSm - 1; x++, index++)
 456           if (obstacles[index] && obstacles[index] != RETIRED)
 457           {
 458             float sum = 0.0f;
 459             int valid = 0;
 460
 461             if (!obstacles[index + 1] || obstacles[index + 1] == RETIRED)
 462             {
 463               sum += _energy[index + 1];
 464               valid++;
 465             }
 466             if (!obstacles[index - 1] || obstacles[index - 1] == RETIRED)
 467             {
 468               sum += _energy[index - 1];
 469               valid++;
 470             }
 471             if (!obstacles[index + _xResSm] || obstacles[index + _xResSm] == RETIRED)
 472             {
 473               sum += _energy[index + _xResSm];
 474               valid++;
 475             }
 476             if (!obstacles[index - _xResSm] || obstacles[index - _xResSm] == RETIRED)
 477             {
 478               sum += _energy[index - _xResSm];
 479               valid++;
 480             }
 481             if (!obstacles[index + _slabSizeSm] || obstacles[index + _slabSizeSm] == RETIRED)
 482             {
 483               sum += _energy[index + _slabSizeSm];
 484               valid++;
 485             }
 486             if (!obstacles[index - _slabSizeSm] || obstacles[index - _slabSizeSm] == RETIRED)
 487             {
 488               sum += _energy[index - _slabSizeSm];
 489               valid++;
 490             }
 491             if (valid > 0)
 492             {
 493               _energy[index] = sum / (float)valid;
 494               obstacles[index] = MARCHED;
 495             }
 496           }
 497     index = _slabSizeSm + _xResSm + 1;
 498     for (int z = 1; z < _zResSm - 1; z++, index += 2 * _xResSm)
 499       for (int y = 1; y < _yResSm - 1; y++, index += 2)
 500         for (int x = 1; x < _xResSm - 1; x++, index++)
 501           if (obstacles[index] == MARCHED)
 502             obstacles[index] = RETIRED;
 503   }
 504   index = _slabSizeSm + _xResSm + 1;
 505   for (int z = 1; z < _zResSm - 1; z++, index += 2 * _xResSm)
 506     for (int y = 1; y < _yResSm - 1; y++, index += 2)
 507       for (int x = 1; x < _xResSm - 1; x++, index++)
 508         if (obstacles[index])
 509           obstacles[index] = 1;
 510 }
 511
 512 //////////////////////////////////////////////////////////////////////////////////////////
 513 // Evaluate derivatives
 514 //////////////////////////////////////////////////////////////////////////////////////////
 515 Vec3 WTURBULENCE::WVelocity(Vec3 orgPos)
 516 {
 517   // arbitrarily offset evaluation points
 518   const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2.0,0,0);
 519   const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2.0,0);
 520   const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2.0);
 521
 522   const float f1y = WNoiseDy(p1, _noiseTile);
 523   const float f1z = WNoiseDz(p1, _noiseTile);
 524
 525   const float f2x = WNoiseDx(p2, _noiseTile);
 526   const float f2z = WNoiseDz(p2, _noiseTile);
 527
 528   const float f3x = WNoiseDx(p3, _noiseTile);
 529   const float f3y = WNoiseDy(p3, _noiseTile);
 530
 531   Vec3 ret = Vec3(
 532       f3y - f2z,
 533       f1z - f3x,
 534       f2x - f1y );
 535   return ret;
 536 }
 537
 538 //////////////////////////////////////////////////////////////////////////////////////////
 539 // Evaluate derivatives with Jacobian
 540 //////////////////////////////////////////////////////////////////////////////////////////
 541 Vec3 WTURBULENCE::WVelocityWithJacobian(Vec3 orgPos, float* xUnwarped, float* yUnwarped, float* zUnwarped)
 542 {
 543   // arbitrarily offset evaluation points
 544   const Vec3 p1 = orgPos + Vec3(NOISE_TILE_SIZE/2.0,0,0);
 545   const Vec3 p2 = orgPos + Vec3(0,NOISE_TILE_SIZE/2.0,0);
 546   const Vec3 p3 = orgPos + Vec3(0,0,NOISE_TILE_SIZE/2.0);
 547
 548   Vec3 final;
 549   final[0] = WNoiseDx(p1, _noiseTile);
 550   final[1] = WNoiseDy(p1, _noiseTile);
 551   final[2] = WNoiseDz(p1, _noiseTile);
 552   // UNUSED const float f1x = xUnwarped[0] * final[0] + xUnwarped[1] * final[1] + xUnwarped[2] * final[2];
 553   const float f1y = yUnwarped[0] * final[0] + yUnwarped[1] * final[1] + yUnwarped[2] * final[2];
 554   const float f1z = zUnwarped[0] * final[0] + zUnwarped[1] * final[1] + zUnwarped[2] * final[2];
 555
 556   final[0] = WNoiseDx(p2, _noiseTile);
 557   final[1] = WNoiseDy(p2, _noiseTile);
 558   final[2] = WNoiseDz(p2, _noiseTile);
 559   const float f2x = xUnwarped[0] * final[0] + xUnwarped[1] * final[1] + xUnwarped[2] * final[2];
 560   // UNUSED const float f2y = yUnwarped[0] * final[0] + yUnwarped[1] * final[1] + yUnwarped[2] * final[2];
 561   const float f2z = zUnwarped[0] * final[0] + zUnwarped[1] * final[1] + zUnwarped[2] * final[2];
 562
 563   final[0] = WNoiseDx(p3, _noiseTile);
 564   final[1] = WNoiseDy(p3, _noiseTile);
 565   final[2] = WNoiseDz(p3, _noiseTile);
 566   const float f3x = xUnwarped[0] * final[0] + xUnwarped[1] * final[1] + xUnwarped[2] * final[2];
 567   const float f3y = yUnwarped[0] * final[0] + yUnwarped[1] * final[1] + yUnwarped[2] * final[2];
 568   // UNUSED const float f3z = zUnwarped[0] * final[0] + zUnwarped[1] * final[1] + zUnwarped[2] * final[2];
 569
 570   Vec3 ret = Vec3(
 571       f3y - f2z,
 572       f1z - f3x,
 573       f2x - f1y );
 574   return ret;
 575 }
 576
 577
 578 //////////////////////////////////////////////////////////////////////
 579 // perform an actual noise advection step
 580 //////////////////////////////////////////////////////////////////////
 581 /*void WTURBULENCE::stepTurbulenceReadable(float dtOrg, float* xvel, float* yvel, float* zvel, unsigned char *obstacles)
 582 {
 583         // enlarge timestep to match grid
 584         const float dt = dtOrg * _amplify;
 585         const float invAmp = 1.0f / _amplify;
 586         float *tempBig1 = new float[_totalCellsBig];
 587         float *tempBig2 = new float[_totalCellsBig];
 588         float *bigUx = new float[_totalCellsBig];
 589         float *bigUy = new float[_totalCellsBig];
 590         float *bigUz = new float[_totalCellsBig];
 591         float *_energy = new float[_totalCellsSm];
 592         float *highFreqEnergy = new float[_totalCellsSm];
 593         float *eigMin  = new float[_totalCellsSm];
 594         float *eigMax  = new float[_totalCellsSm];
 595
 596         memset(tempBig1, 0, sizeof(float)*_totalCellsBig);
 597         memset(tempBig2, 0, sizeof(float)*_totalCellsBig);
 598         memset(highFreqEnergy, 0, sizeof(float)*_totalCellsSm);
 599         memset(eigMin, 0, sizeof(float)*_totalCellsSm);
 600         memset(eigMax, 0, sizeof(float)*_totalCellsSm);
 601
 602         // prepare textures
 603         advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempBig1, tempBig2);
 604
 605   // compute eigenvalues of the texture coordinates
 606   computeEigenvalues(eigMin, eigMax);
 607
 608   // do wavelet decomposition of energy
 609   computeEnergy(_energy, xvel, yvel, zvel, obstacles);
 610   decomposeEnergy(_energy, highFreqEnergy);
 611
 612   // zero out coefficients inside of the obstacle
 613   for (int x = 0; x < _totalCellsSm; x++)
 614     if (obstacles[x]) _energy[x] = 0.f;
 615
 616   float maxVelocity = 0.;
 617   for (int z = 1; z < _zResBig - 1; z++)
 618     for (int y = 1; y < _yResBig - 1; y++)
 619       for (int x = 1; x < _xResBig - 1; x++)
 620       {
 621         // get unit position for both fine and coarse grid
 622         const Vec3 pos = Vec3(x,y,z);
 623         const Vec3 posSm = pos * invAmp;
 624
 625         // get grid index for both fine and coarse grid
 626         const int index = x + y *_xResBig + z *_slabSizeBig;
 627         const int indexSmall = (int)posSm[0] + (int)posSm[1] * _xResSm + (int)posSm[2] * _slabSizeSm;
 628
 629         // get a linearly interpolated velocity and texcoords
 630         // from the coarse grid
 631         Vec3 vel = INTERPOLATE::lerp3dVec( xvel,yvel,zvel,
 632             posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
 633         Vec3 uvw = INTERPOLATE::lerp3dVec( _tcU,_tcV,_tcW,
 634             posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
 635
 636         // multiply the texture coordinate by _resSm so that turbulence
 637         // synthesis begins at the first octave that the coarse grid
 638         // cannot capture
 639         Vec3 texCoord = Vec3(uvw[0] * _resSm[0],
 640                              uvw[1] * _resSm[1],
 641                              uvw[2] * _resSm[2]);
 642
 643         // retrieve wavelet energy at highest frequency
 644         float energy = INTERPOLATE::lerp3d(
 645             highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm);
 646
 647         // base amplitude for octave 0
 648         float coefficient = sqrtf(2.0f * fabs(energy));
 649         const float amplitude = *_strength * fabs(0.5 * coefficient) * persistence;
 650
 651         // add noise to velocity, but only if the turbulence is
 652         // sufficiently undeformed, and the energy is large enough
 653         // to make a difference
 654         const bool addNoise = eigMax[indexSmall] < 2. &&
 655                               eigMin[indexSmall] > 0.5;
 656         if (addNoise && amplitude > _cullingThreshold) {
 657           // base amplitude for octave 0
 658           float amplitudeScaled = amplitude;
 659
 660           for (int octave = 0; octave < _octaves; octave++)
 661           {
 662             // multiply the vector noise times the maximum allowed
 663             // noise amplitude at this octave, and add it to the total
 664             vel += WVelocity(texCoord) * amplitudeScaled;
 665
 666             // scale coefficient for next octave
 667             amplitudeScaled *= persistence;
 668             texCoord *= 2.0f;
 669           }
 670         }
 671
 672         // Store velocity + turbulence in big grid for maccormack step
 673         //
 674         // If you wanted to save memory, you would instead perform a
 675         // semi-Lagrangian backtrace for the current grid cell here. Then
 676         // you could just throw the velocity away.
 677         bigUx[index] = vel[0];
 678         bigUy[index] = vel[1];
 679         bigUz[index] = vel[2];
 680
 681         // compute the velocity magnitude for substepping later
 682         const float velMag = bigUx[index] * bigUx[index] +
 683                              bigUy[index] * bigUy[index] +
 684                              bigUz[index] * bigUz[index];
 685         if (velMag > maxVelocity) maxVelocity = velMag;
 686
 687         // zero out velocity inside obstacles
 688         float obsCheck = INTERPOLATE::lerp3dToFloat(
 689             obstacles, posSm[0], posSm[1], posSm[2], _xResSm, _yResSm, _zResSm);
 690         if (obsCheck > 0.95)
 691           bigUx[index] = bigUy[index] = bigUz[index] = 0.;
 692       }
 693
 694   // prepare density for an advection
 695   SWAP_POINTERS(_densityBig, _densityBigOld);
 696
 697   // based on the maximum velocity present, see if we need to substep,
 698   // but cap the maximum number of substeps to 5
 699   const int maxSubSteps = 5;
 700   maxVelocity = sqrt(maxVelocity) * dt;
 701   int totalSubsteps = (int)(maxVelocity / (float)maxSubSteps);
 702   totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps;
 703   totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps;
 704   const float dtSubdiv = dt / (float)totalSubsteps;
 705
 706   // set boundaries of big velocity grid
 707   FLUID_3D::setZeroX(bigUx, _resBig, 0, _resBig[2]);
 708   FLUID_3D::setZeroY(bigUy, _resBig, 0, _resBig[2]);
 709   FLUID_3D::setZeroZ(bigUz, _resBig, 0, _resBig[2]);
 710
 711   // do the MacCormack advection, with substepping if necessary
 712   for(int substep = 0; substep < totalSubsteps; substep++)
 713   {
 714     FLUID_3D::advectFieldMacCormack(dtSubdiv, bigUx, bigUy, bigUz,
 715         _densityBigOld, _densityBig, tempBig1, tempBig2, _resBig, NULL);
 716
 717     if (substep < totalSubsteps - 1)
 718       SWAP_POINTERS(_densityBig, _densityBigOld);
 719   } // substep
 720
 721   // wipe the density borders
 722   FLUID_3D::setZeroBorder(_densityBig, _resBig, 0, _resBig[2]);
 723
 724   // reset texture coordinates now in preparation for next timestep
 725   // Shouldn't do this before generating the noise because then the
 726   // eigenvalues stored do not reflect the underlying texture coordinates
 727   resetTextureCoordinates(eigMin, eigMax);
 728
 729   delete[] tempBig1;
 730   delete[] tempBig2;
 731   delete[] bigUx;
 732   delete[] bigUy;
 733   delete[] bigUz;
 734   delete[] _energy;
 735   delete[] highFreqEnergy;
 736
 737   delete[] eigMin;
 738   delete[] eigMax;
 739
 740
 741   _totalStepsBig++;
 742 }*/
 743
 744 //struct
 745
 746 //////////////////////////////////////////////////////////////////////
 747 // perform the full turbulence algorithm, including OpenMP
 748 // if available
 749 //////////////////////////////////////////////////////////////////////
 750 void WTURBULENCE::stepTurbulenceFull(float dtOrg, float* xvel, float* yvel, float* zvel, unsigned char *obstacles)
 751 {
 752         // enlarge timestep to match grid
 753         const float dt = dtOrg * _amplify;
 754         const float invAmp = 1.0f / _amplify;
 755         float *tempBig1 = (float *)calloc(_totalCellsBig, sizeof(float));
 756         float *tempBig2 = (float *)calloc(_totalCellsBig, sizeof(float));
 757         float *bigUx = (float *)calloc(_totalCellsBig, sizeof(float));
 758         float *bigUy = (float *)calloc(_totalCellsBig, sizeof(float));
 759         float *bigUz = (float *)calloc(_totalCellsBig, sizeof(float));
 760         float *_energy = (float *)calloc(_totalCellsSm, sizeof(float));
 761         float *highFreqEnergy = (float *)calloc(_totalCellsSm, sizeof(float));
 762         float *eigMin  = (float *)calloc(_totalCellsSm, sizeof(float));
 763         float *eigMax  = (float *)calloc(_totalCellsSm, sizeof(float));
 764
 765         memset(_tcTemp, 0, sizeof(float)*_totalCellsSm);
 766
 767
 768         // prepare textures
 769         advectTextureCoordinates(dtOrg, xvel,yvel,zvel, tempBig1, tempBig2);
 770
 771         // do wavelet decomposition of energy
 772         computeEnergy(_energy, xvel, yvel, zvel, obstacles);
 773
 774         for (int x = 0; x < _totalCellsSm; x++)
 775                 if (obstacles[x]) _energy[x] = 0.f;
 776
 777         decomposeEnergy(_energy, highFreqEnergy);
 778
 779         // zero out coefficients inside of the obstacle
 780         for (int x = 0; x < _totalCellsSm; x++)
 781                 if (obstacles[x]) highFreqEnergy[x] = 0.f;
 782
 783         Vec3Int ressm(_xResSm, _yResSm, _zResSm);
 784         FLUID_3D::setNeumannX(highFreqEnergy, ressm, 0 , ressm[2]);
 785         FLUID_3D::setNeumannY(highFreqEnergy, ressm, 0 , ressm[2]);
 786         FLUID_3D::setNeumannZ(highFreqEnergy, ressm, 0 , ressm[2]);
 787
 788
 789    int threadval = 1;
 790 #if PARALLEL==1
 791   threadval = omp_get_max_threads();
 792 #endif
 793
 794
 795   // parallel region setup
 796   // Uses omp_get_max_trheads to get number of required cells.
 797   float* maxVelMagThreads = new float[threadval];
 798
 799   for (int i=0; i<threadval; i++) maxVelMagThreads[i] = -1.0f;
 800
 801 #if PARALLEL==1
 802
 803 #pragma omp parallel
 804 #endif
 805   { float maxVelMag1 = 0.;
 806 #if PARALLEL==1
 807     const int id  = omp_get_thread_num(); /*, num = omp_get_num_threads(); */
 808 #endif
 809
 810   // vector noise main loop
 811 #if PARALLEL==1
 812 #pragma omp for schedule(static,1)
 813 #endif
 814   for (int zSmall = 0; zSmall < _zResSm; zSmall++)
 815   {
 816   for (int ySmall = 0; ySmall < _yResSm; ySmall++)
 817   for (int xSmall = 0; xSmall < _xResSm; xSmall++)
 818   {
 819     const int indexSmall = xSmall + ySmall * _xResSm + zSmall * _slabSizeSm;
 820
 821     // compute jacobian
 822     float jacobian[3][3] = {
 823       { minDx(xSmall, ySmall, zSmall, _tcU, _resSm), minDx(xSmall, ySmall, zSmall, _tcV, _resSm), minDx(xSmall, ySmall, zSmall, _tcW, _resSm) } ,
 824       { minDy(xSmall, ySmall, zSmall, _tcU, _resSm), minDy(xSmall, ySmall, zSmall, _tcV, _resSm), minDy(xSmall, ySmall, zSmall, _tcW, _resSm) } ,
 825       { minDz(xSmall, ySmall, zSmall, _tcU, _resSm), minDz(xSmall, ySmall, zSmall, _tcV, _resSm), minDz(xSmall, ySmall, zSmall, _tcW, _resSm) }
 826     };
 827
 828     // get LU factorization of texture jacobian and apply
 829     // it to unit vectors
 830     sLU LU = computeLU(jacobian);
 831     float xUnwarped[3], yUnwarped[3], zUnwarped[3];
 832     float xWarped[3], yWarped[3], zWarped[3];
 833     bool nonSingular = isNonsingular(LU);
 834
 835         xUnwarped[0] = 1.0f; xUnwarped[1] = 0.0f; xUnwarped[2] = 0.0f;
 836         yUnwarped[0] = 0.0f; yUnwarped[1] = 1.0f; yUnwarped[2] = 0.0f;
 837         zUnwarped[0] = 0.0f; zUnwarped[1] = 0.0f; zUnwarped[2] = 1.0f;
 838
 839         xWarped[0] = 1.0f; xWarped[1] = 0.0f; xWarped[2] = 0.0f;
 840         yWarped[0] = 0.0f; yWarped[1] = 1.0f; yWarped[2] = 0.0f;
 841         zWarped[0] = 0.0f; zWarped[1] = 0.0f; zWarped[2] = 1.0f;
 842
 843 #if 0
 844         // UNUSED
 845     float eigMax = 10.0f;
 846     float eigMin = 0.1f;
 847 #endif
 848     if (nonSingular)
 849     {
 850       solveLU3x3(LU, xUnwarped, xWarped);
 851       solveLU3x3(LU, yUnwarped, yWarped);
 852       solveLU3x3(LU, zUnwarped, zWarped);
 853
 854       // compute the eigenvalues while we have the Jacobian available
 855       Vec3 eigenvalues = Vec3(1.);
 856       computeEigenvalues3x3( &eigenvalues[0], jacobian);
 857       eigMax[indexSmall] = MAX3V(eigenvalues);
 858       eigMin[indexSmall] = MIN3V(eigenvalues);
 859     }
 860
 861     // make sure to skip one on the beginning and end
 862     int xStart = (xSmall == 0) ? 1 : 0;
 863     int xEnd   = (xSmall == _xResSm - 1) ? _amplify - 1 : _amplify;
 864     int yStart = (ySmall == 0) ? 1 : 0;
 865     int yEnd   = (ySmall == _yResSm - 1) ? _amplify - 1 : _amplify;
 866     int zStart = (zSmall == 0) ? 1 : 0;
 867     int zEnd   = (zSmall == _zResSm - 1) ? _amplify - 1 : _amplify;
 868
 869     for (int zBig = zStart; zBig < zEnd; zBig++)
 870     for (int yBig = yStart; yBig < yEnd; yBig++)
 871     for (int xBig = xStart; xBig < xEnd; xBig++)
 872     {
 873       const int x = xSmall * _amplify + xBig;
 874       const int y = ySmall * _amplify + yBig;
 875       const int z = zSmall * _amplify + zBig;
 876
 877       // get unit position for both fine and coarse grid
 878       const Vec3 pos = Vec3(x,y,z);
 879       const Vec3 posSm = pos * invAmp;
 880
 881       // get grid index for both fine and coarse grid
 882       const int index = x + y *_xResBig + z *_slabSizeBig;
 883
 884       // get a linearly interpolated velocity and texcoords
 885       // from the coarse grid
 886       Vec3 vel = INTERPOLATE::lerp3dVec( xvel,yvel,zvel,
 887           posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
 888       Vec3 uvw = INTERPOLATE::lerp3dVec( _tcU,_tcV,_tcW,
 889           posSm[0], posSm[1], posSm[2], _xResSm,_yResSm,_zResSm);
 890
 891       // multiply the texture coordinate by _resSm so that turbulence
 892       // synthesis begins at the first octave that the coarse grid
 893       // cannot capture
 894       Vec3 texCoord = Vec3(uvw[0] * _resSm[0],
 895                            uvw[1] * _resSm[1],
 896                            uvw[2] * _resSm[2]);
 897
 898       // retrieve wavelet energy at highest frequency
 899       float energy = INTERPOLATE::lerp3d(
 900           highFreqEnergy, posSm[0],posSm[1],posSm[2], _xResSm, _yResSm, _zResSm);
 901
 902       // base amplitude for octave 0
 903       float coefficient = sqrtf(2.0f * fabs(energy));
 904       const float amplitude = *_strength * fabs(0.5 * coefficient) * persistence;
 905
 906       // add noise to velocity, but only if the turbulence is
 907       // sufficiently undeformed, and the energy is large enough
 908       // to make a difference
 909       const bool addNoise = eigMax[indexSmall] < 2. &&
 910                             eigMin[indexSmall] > 0.5;
 911       if (addNoise && amplitude > _cullingThreshold) {
 912         // base amplitude for octave 0
 913         float amplitudeScaled = amplitude;
 914
 915         for (int octave = 0; octave < _octaves; octave++)
 916         {
 917           // multiply the vector noise times the maximum allowed
 918           // noise amplitude at this octave, and add it to the total
 919           vel += WVelocityWithJacobian(texCoord, &xUnwarped[0], &yUnwarped[0], &zUnwarped[0]) * amplitudeScaled;
 920
 921           // scale coefficient for next octave
 922           amplitudeScaled *= persistence;
 923           texCoord *= 2.0f;
 924         }
 925       }
 926
 927       // Store velocity + turbulence in big grid for maccormack step
 928       //
 929       // If you wanted to save memory, you would instead perform a
 930       // semi-Lagrangian backtrace for the current grid cell here. Then
 931       // you could just throw the velocity away.
 932       bigUx[index] = vel[0];
 933       bigUy[index] = vel[1];
 934       bigUz[index] = vel[2];
 935
 936       // compute the velocity magnitude for substepping later
 937       const float velMag = bigUx[index] * bigUx[index] +
 938                            bigUy[index] * bigUy[index] +
 939                            bigUz[index] * bigUz[index];
 940       if (velMag > maxVelMag1) maxVelMag1 = velMag;
 941
 942       // zero out velocity inside obstacles
 943       float obsCheck = INTERPOLATE::lerp3dToFloat(
 944           obstacles, posSm[0], posSm[1], posSm[2], _xResSm, _yResSm, _zResSm);
 945       if (obsCheck > 0.95)
 946         bigUx[index] = bigUy[index] = bigUz[index] = 0.;
 947     } // xyz*/
 948
 949 #if PARALLEL==1
 950     maxVelMagThreads[id] = maxVelMag1;
 951 #else
 952     maxVelMagThreads[0] = maxVelMag1;
 953 #endif
 954   }
 955   }
 956   } // omp
 957
 958   // compute maximum over threads
 959   float maxVelMag = maxVelMagThreads[0];
 960 #if PARALLEL==1
 961   for (int i = 1; i < threadval; i++)
 962     if (maxVelMag < maxVelMagThreads[i])
 963       maxVelMag = maxVelMagThreads[i];
 964 #endif
 965   delete [] maxVelMagThreads;
 966
 967
 968   // prepare density for an advection
 969   SWAP_POINTERS(_densityBig, _densityBigOld);
 970
 971   // based on the maximum velocity present, see if we need to substep,
 972   // but cap the maximum number of substeps to 5
 973   const int maxSubSteps = 25;
 974   const int maxVel = 5;
 975   maxVelMag = sqrt(maxVelMag) * dt;
 976   int totalSubsteps = (int)(maxVelMag / (float)maxVel);
 977   totalSubsteps = (totalSubsteps < 1) ? 1 : totalSubsteps;
 978   // printf("totalSubsteps: %d\n", totalSubsteps);
 979   totalSubsteps = (totalSubsteps > maxSubSteps) ? maxSubSteps : totalSubsteps;
 980   const float dtSubdiv = dt / (float)totalSubsteps;
 981
 982   // set boundaries of big velocity grid
 983   FLUID_3D::setZeroX(bigUx, _resBig, 0 , _resBig[2]);
 984   FLUID_3D::setZeroY(bigUy, _resBig, 0 , _resBig[2]);
 985   FLUID_3D::setZeroZ(bigUz, _resBig, 0 , _resBig[2]);
 986
 987 #if PARALLEL==1
 988   int stepParts = threadval*2;  // Dividing parallelized sections into numOfThreads * 2 sections
 989   float partSize = (float)_zResBig/stepParts;   // Size of one part;
 990
 991   if (partSize < 4) {stepParts = threadval;                                     // If the slice gets too low (might actually slow things down, change it to larger
 992                                         partSize = (float)_zResBig/stepParts;}
 993   if (partSize < 4) {stepParts = (int)(ceil((float)_zResBig/4.0f));     // If it's still too low (only possible on future systems with +24 cores), change it to 4
 994                                         partSize = (float)_zResBig/stepParts;}
 995 #else
 996   int zBegin=0;
 997   int zEnd=_resBig[2];
 998 #endif
 999
1000   // do the MacCormack advection, with substepping if necessary
1001   for(int substep = 0; substep < totalSubsteps; substep++)
1002   {
1003
1004 #if PARALLEL==1
1005         #pragma omp parallel
1006         {
1007
1008         #pragma omp for schedule(static,1)
1009         for (int i=0; i<stepParts; i++)
1010         {
1011                 int zBegin = (int)((float)i*partSize + 0.5f);
1012                 int zEnd = (int)((float)(i+1)*partSize + 0.5f);
1013 #endif
1014                 FLUID_3D::advectFieldMacCormack1(dtSubdiv, bigUx, bigUy, bigUz,
1015                     _densityBigOld, tempBig1, _resBig, zBegin, zEnd);
1016 #if PARALLEL==1
1017         }
1018
1019         #pragma omp barrier
1020
1021         #pragma omp for schedule(static,1)
1022         for (int i=0; i<stepParts; i++)
1023         {
1024                 int zBegin = (int)((float)i*partSize + 0.5f);
1025                 int zEnd = (int)((float)(i+1)*partSize + 0.5f);
1026 #endif
1027                 FLUID_3D::advectFieldMacCormack2(dtSubdiv, bigUx, bigUy, bigUz,
1028                     _densityBigOld, _densityBig, tempBig1, tempBig2, _resBig, NULL, zBegin, zEnd);
1029 #if PARALLEL==1
1030         }
1031         }
1032 #endif
1033
1034     if (substep < totalSubsteps - 1)
1035       SWAP_POINTERS(_densityBig, _densityBigOld);
1036   } // substep
1037
1038   free(tempBig1);
1039   free(tempBig2);
1040   free(bigUx);
1041   free(bigUy);
1042   free(bigUz);
1043   free(_energy);
1044   free(highFreqEnergy);
1045
1046   // wipe the density borders
1047   FLUID_3D::setZeroBorder(_densityBig, _resBig, 0 , _resBig[2]);
1048
1049   // reset texture coordinates now in preparation for next timestep
1050   // Shouldn't do this before generating the noise because then the
1051   // eigenvalues stored do not reflect the underlying texture coordinates
1052   resetTextureCoordinates(eigMin, eigMax);
1053
1054   free(eigMin);
1055   free(eigMax);
1056
1057   // output files
1058   // string prefix = string("./amplified.preview/density_bigxy_");
1059   // FLUID_3D::writeImageSliceXY(_densityBig, _resBig, _resBig[2]/2, prefix, _totalStepsBig, 1.0f);
1060   //string df3prefix = string("./df3/density_big_");
1061   //IMAGE::dumpDF3(_totalStepsBig, df3prefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]);
1062   // string pbrtPrefix = string("./pbrt/density_big_");
1063   // IMAGE::dumpPBRT(_totalStepsBig, pbrtPrefix, _densityBig, _resBig[0],_resBig[1],_resBig[2]);
1064
1065   _totalStepsBig++;
1066 }