1 /* -*- c++ -*- */ 2 #ifndef MULTIGRID_H 3 #define MULTIGRID_H 4 5 #include "Array.h" 6 #include "Vector3D.h" 7 #include "mathutilities.h" 8 #include "Report.h" 9 10 //#define DEBUG_MULTIGRID 11 //#define DEBUG_MULTIGRID_TIMING 12 13 namespace ProtoMol { 14 //_________________________________________________________________ MultiGrid 15 static const Real BORDER_TOLERANCE = 0.0001; 16 static const int BORDER = 2; 17 18 19 /** 20 * Multi grid algorithm based on TInterpolation scheme using TKernel 21 * handling non-periodic and periodic boundary conditions. 22 */ 23 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> 24 class MultiGrid { 25 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 26 // Typedef & const 27 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 28 private: 29 struct Int3D {int x, y, z;}; 30 /// 3D interpolation 31 struct Interpolation3D { 32 TInterpolation x; 33 TInterpolation y; 34 TInterpolation z; Interpolation3DInterpolation3D35 Interpolation3D(){}; Interpolation3DInterpolation3D36 Interpolation3D(unsigned int order):x(order),y(order),z(order){}; 37 }; 38 39 40 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 41 // Constructors, destructors, assignment 42 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 43 public: 44 MultiGrid(); 45 #ifdef DEBUG_MULTIGRID_TIMING ~MultiGrid()46 ~MultiGrid(){ 47 if(myCounterDirect+myCounterCorrection > 0) 48 Report::report << allnodes << plain 49 << "[MultiGrid] Smooth part: direct=" 50 << myCounterDirect<<", correction=" 51 << myCounterCorrection<<"."<<Report::endr; 52 } getCounter()53 long getCounter(){return (myCounterDirect + myCounterCorrection);} 54 #endif 55 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 56 // New methods of class MultiGrid 57 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 58 public: 59 60 /// Interpolation weights for the position and anterpolate the charge q onto grid level 0, q -> Q(0) 61 void anterpolateCharge(Real q, const Vector3D& coord, unsigned int index); 62 63 /// Anterpolates Q from level to level+1, Q(level) -> Q(level+1) 64 void fineToCoarse(int level, unsigned int block, unsigned int nBlocks); fineToCoarse(int level)65 void fineToCoarse(int level){fineToCoarse(level,0,1);} 66 67 /// Computes the potential for the top level, Q(maxLevels-1) -> V(maxLevels-1) 68 void direct(unsigned int block, unsigned int nBlocks); direct()69 void direct(){direct(0,1);}; 70 71 /// Interpolates Q from level to level-1, Q(level) -> Q(level-1) 72 void coarseToFine(int level, unsigned int block, unsigned int nBlocks); coarseToFine(int level)73 void coarseToFine(int level){coarseToFine(level,0,1);} 74 75 /// Adds the correction term Q(level) -> V(level) 76 void correction(int level, unsigned int block, unsigned int nBlocks); correction(int level)77 void correction(int level){correction(level,0,1);}; 78 79 /// Computes the energy of a given level, V(level) -> energy 80 Real energy(int level, unsigned int block, unsigned int nBlocks); energy(int level)81 Real energy(int level){return energy(level,0,1);}; 82 83 /// Interpolates the force from grid, level 0, Q(0) -> force 84 void interpolateForce(Real q, unsigned int index, Vector3D& force); 85 86 87 void initialize(unsigned int n, Real s, int levels, 88 int nx, int ny, int nz, 89 int interOrder, int ratio, 90 Vector3D min, Vector3D max, 91 Vector3D h, Vector3D origin); 92 93 /** Update the length, width and height of the grid according min and max such 94 * that anter-/interpolated works regardless order and boundary conditions. 95 * In case of non-PBC it adds a border such that all points inside [min,max] 96 * are correctly anter-/interpolated. 97 */ 98 void updateSize(Vector3D min, Vector3D max); 99 getV(int level,Real * & begin,Real * & end)100 void getV(int level, Real*& begin, Real*& end) {begin=myV[level].begin();end=myV[level].end();} getQ(int level,Real * & begin,Real * & end)101 void getQ(int level, Real*& begin, Real*& end) {begin=myQ[level].begin();end=myQ[level].end();} 102 void clear(); 103 104 Real sumV(int level); 105 Real sumQ(int level); print(int level)106 void print(int level){printQ(level);printV(level);} 107 void printV(int level); 108 void printQ(int level); 109 void printConst(); 110 private: 111 void precomputeG(); 112 void blocksCorrection(unsigned int nBlocks); 113 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 114 // My data members 115 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 116 private: 117 std::vector<Array<Real,3> > myQ; ///< Arrays of charges 118 std::vector<Array<Real,3> > myV; ///< Arrays of potentials 119 Array<Real,3> myGDirect; ///< G_direct for top level, G->V 120 Array<Real,3> myGCorrection; ///< G_Correction for other levels, G->V 121 int myGCorrDimX; ///< Dimensions of G_correction 122 int myGCorrDimY; 123 int myGCorrDimZ; 124 Vector3D myMin; // 125 Vector3D myMax; // 126 Vector3D myCMin; // 127 Vector3D myCMax; // 128 Vector3D myH; // 129 Vector3D myOrigin; // 130 int myNX; ///< x dimension of the grid level 0 131 int myNY; ///< y dimension of the grid level 0 132 int myNZ; ///< z dimension of the grid level 0 133 int myNXOffset; ///< x index offset between to grid levels 134 int myNYOffset; ///< y index offset between to grid levels 135 int myNZOffset; ///< z index offset between to grid levels 136 Real myHX; ///< h_x at level 0 137 Real myHY; ///< h_y at level 0 138 Real myHZ; ///< h_z at level 0 139 Real myHXr; ///< 1/h_x at level 0 140 Real myHYr; ///< 1/h_y at level 0 141 Real myHZr; ///< 1/h_z at level 0 142 std::vector<Int3D> myScaledParticleIntPositions; ///< Starting index for anter/-interpolation 143 // for each particle 144 std::vector<Interpolation3D> myInterpolations; ///< Interpolation weights for each particle 145 std::vector<TInterpolation> myGridInterpolation; ///< Interpolation between grids 146 Real myS; ///< Softening distance 147 Real myRS; ///< 1 / Softening distance 148 int myInterOrder; ///< Interpolation order 149 int myRatio; ///< Ratio between to levels, usually 2 150 int myLevels; ///< Number of levels, 0,1,2,...,myLevels-1 151 int mySigma; ///< 1 if interpolation weights ...,0 for w=0, 0 else 152 std::vector<Int3D> myDim; ///< Dimension of the grids 153 std::vector<Real> myScale; ///< Scaling factor for each level 154 std::vector<std::vector<int> > myBlocksCorrection; ///< List of blocks for the correction 155 #ifdef DEBUG_MULTIGRID_TIMING 156 long myCounterDirect; 157 long myCounterCorrection; 158 #endif 159 }; 160 //______________________________________________________________________ INLINES 161 162 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> anterpolateCharge(Real q,const Vector3D & coord,unsigned int index)163 inline void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::anterpolateCharge(Real q, const Vector3D& coord, unsigned int index) { 164 Real x = (coord.x-myMin.x)*myHXr; 165 Real y = (coord.y-myMin.y)*myHYr; 166 Real z = (coord.z-myMin.z)*myHZr; 167 if(pbcX) while(x < 0.0) x += myNX; 168 if(pbcX) while(x >= myNX) x -= myNX; 169 if(pbcY) while(y < 0.0) y += myNY; 170 if(pbcY) while(y >= myNY) y -= myNY; 171 if(pbcZ) while(z < 0.0) z += myNZ; 172 if(pbcZ) while(z >= myNZ) z -= myNZ; 173 int intX = (int)x; 174 int intY = (int)y; 175 int intZ = (int)z; 176 int i0 = intX+myNXOffset; 177 int j0 = intY+myNYOffset; 178 int k0 = intZ+myNZOffset; 179 myScaledParticleIntPositions[index].x = i0; 180 myScaledParticleIntPositions[index].y = j0; 181 myScaledParticleIntPositions[index].z = k0; 182 myInterpolations[index].x.set(x-intX); 183 myInterpolations[index].y.set(y-intY); 184 myInterpolations[index].z.set(z-intZ); 185 186 Real* thetaX = myInterpolations[index].x.theta; 187 Real* thetaY = myInterpolations[index].y.theta; 188 Real* thetaZ = myInterpolations[index].z.theta; 189 Array<Real,3>& rQ = myQ[0]; 190 for(int i=0;i<myInterOrder;i++){ 191 Real a = q*thetaX[i]; 192 int i1 = i+i0; 193 if(pbcX) i1 = i1 % myNX; 194 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 195 Array<Real,3>::RefArray<2> rQX = rQ[i1]; 196 #else 197 RefArray<Real,2> rQX = rQ[i1]; 198 #endif 199 for(int j=0;j<myInterOrder;j++){ 200 Real ab = a*thetaY[j]; 201 int j1 = j+j0; 202 if(pbcY) j1 = j1 % myNY; 203 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 204 Array<Real,3>::RefArray<1> rQXY = rQX[j1]; 205 #else 206 RefArray<Real,1> rQXY = rQX[j1]; 207 #endif 208 for(int k=0;k<myInterOrder;k++){ 209 int k1 = k+k0; 210 if(pbcZ) k1 = k1 % myNZ; 211 rQXY[k1] += ab*thetaZ[k]; 212 } 213 } 214 } 215 } 216 217 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> fineToCoarse(int level,unsigned int block,unsigned int nBlocks)218 inline void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::fineToCoarse(int level, unsigned int block, unsigned int nBlocks){ 219 int nx = myDim[level].x; 220 int ny = myDim[level].y; 221 int nz = myDim[level].z; 222 int nx2 = myDim[level+1].x; 223 int ny2 = myDim[level+1].y; 224 int nz2 = myDim[level+1].z; 225 Array<Real,3>& rQ0 = myQ[level]; 226 Array<Real,3>& rQ1 = myQ[level+1]; 227 228 int nyz = ny*nz; 229 int n = nx*nyz; 230 int sn = (n*block)/nBlocks; 231 int en = (n*(block+1))/nBlocks - 1; 232 int count = 0; 233 int size = en-sn+1; 234 if(size == 0) 235 return; 236 237 int k = sn % nz; 238 int j = (sn / nz) % ny; 239 int i = (sn / nyz); 240 int ez = (en % nz)+1; 241 int ey = ((en / nz) % ny)+1; 242 int ex = (en / nyz)+1; 243 if(j < ey-1) 244 ez = nz; 245 if(i < ex-1){ 246 ey = ny; 247 ez = nz; 248 } 249 250 for (; i < ex; i++,j=0){ 251 int i1 = i/myRatio; 252 Real* thetaX = myGridInterpolation[i%myRatio].theta; 253 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 254 Array<Real,3>::RefArray<2> rQ0X = rQ0[i]; 255 #else 256 RefArray<Real,2> rQ0X = rQ0[i]; 257 #endif 258 for (; j < ey; j++,k=0){ 259 int j1 = j/myRatio; 260 Real* thetaY = myGridInterpolation[j%myRatio].theta; 261 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 262 Array<Real,3>::RefArray<1> rQ0XY = rQ0X[j]; 263 #else 264 RefArray<Real,1> rQ0XY = rQ0X[j]; 265 #endif 266 for (; k < ez; k++){ 267 int k1 = k/myRatio; 268 Real* thetaZ = myGridInterpolation[k%myRatio].theta; 269 Real q = rQ0XY[k]; 270 for(int i0=0;i0<myInterOrder;i0++){ 271 int i2 = i1+i0; 272 if(pbcX) i2 = (myNXOffset+i2) % nx2; 273 Real qx = q*thetaX[i0]; 274 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 275 Array<Real,3>::RefArray<2> rQ1X = rQ1[i2]; 276 #else 277 RefArray<Real,2> rQ1X = rQ1[i2]; 278 #endif 279 for(int j0=0;j0<myInterOrder;j0++){ 280 int j2 = j1+j0; 281 if(pbcY) j2 = (myNYOffset+j2) % ny2; 282 Real qxy = qx*thetaY[j0]; 283 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 284 Array<Real,3>::RefArray<1> rQ1XY = rQ1X[j2]; 285 #else 286 RefArray<Real,1> rQ1XY = rQ1X[j2]; 287 #endif 288 for(int k0=0;k0<myInterOrder;k0++){ 289 int k2 = k1+k0; 290 if(pbcZ) k2 = (myNZOffset+k2) % nz2; 291 rQ1XY[k2] += qxy*thetaZ[k0]; 292 } 293 } 294 } 295 count++; 296 if(count >= size) 297 return; 298 } 299 } 300 } 301 } 302 303 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> direct(unsigned int block,unsigned int nBlocks)304 inline void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::direct(unsigned int block, unsigned int nBlocks){ 305 int nx = myDim[myLevels-1].x; 306 int ny = myDim[myLevels-1].y; 307 int nz = myDim[myLevels-1].z; 308 Array<Real,3>& rV = myV[myLevels-1]; 309 Array<Real,3>& rQ = myQ[myLevels-1]; 310 Real g0 = myGDirect[0][0][0]; 311 int nyz = ny*nz; 312 int n = nx*nyz; 313 int sn = 0; 314 int en = n -1; 315 if(block > 0) 316 sn = (int)((-1.0+2.0*n-sqrt(power<2>(1.0-2.0*n)-4.0*n*(n-1.0)*(block)/(Real)nBlocks))/2.0+0.5); 317 if(block < nBlocks-1) 318 en = (int)((-1.0+2.0*n-sqrt(power<2>(1.0-2.0*n)-4.0*n*(n-1.0)*(block+1.0)/(Real)nBlocks))/2.0 - 0.5); 319 320 int count = 0; 321 int size = en-sn+1; 322 if(size == 0) 323 return; 324 325 int k = sn % nz; 326 int j = (sn / nz) % ny; 327 int i = (sn / nyz); 328 int ez = (en % nz)+1; 329 int ey = ((en / nz) % ny)+1; 330 int ex = (en / nyz)+1; 331 if(j < ey-1) 332 ez = nz; 333 if(i < ex-1){ 334 ey = ny; 335 ez = nz; 336 } 337 338 for (; i < ex; i++,j=0){ 339 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 340 Array<Real,3>::RefArray<2> rV1X = rV[i]; 341 Array<Real,3>::RefArray<2> rQ1X = rQ[i]; 342 #else 343 RefArray<Real,2> rV1X = rV[i]; 344 RefArray<Real,2> rQ1X = rQ[i]; 345 #endif 346 for (; j < ey; j++,k=0){ 347 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 348 Array<Real,3>::RefArray<1> rV1XY = rV1X[j]; 349 Array<Real,3>::RefArray<1> rQ1XY = rQ1X[j]; 350 #else 351 RefArray<Real,1> rV1XY = rV1X[j]; 352 RefArray<Real,1> rQ1XY = rQ1X[j]; 353 #endif 354 for (; k < ez; k++){ 355 int l,m,n; 356 Real v = 0.0; 357 Real q = rQ1XY[k]; 358 for (l = i, m = j, n = k+1; l < nx; l++, m=0){ 359 int i0 = i-l; 360 if(pbcX){ 361 i0 = (i0+nx)%nx; 362 i0 = std::min(i0,nx-i0); 363 } 364 else { 365 i0 = i0 < 0 ? -i0 : i0; 366 } 367 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 368 Array<Real,3>::RefArray<2> rGDX = myGDirect[i0]; 369 Array<Real,3>::RefArray<2> rVX = rV[l]; 370 Array<Real,3>::RefArray<2> rQX = rQ[l]; 371 #else 372 RefArray<Real,2> rGDX = myGDirect[i0]; 373 RefArray<Real,2> rVX = rV[l]; 374 RefArray<Real,2> rQX = rQ[l]; 375 #endif 376 for (; m < ny; m++,n=0){ 377 int j0 = j-m; 378 if(pbcX){ 379 j0 = (j0+ny)%ny; 380 j0 = std::min(j0,ny-j0); 381 } 382 else { 383 j0 = j0 < 0 ? -j0 : j0; 384 } 385 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 386 Array<Real,3>::RefArray<1> rGDXY = rGDX[j0]; 387 Array<Real,3>::RefArray<1> rVXY = rVX[m]; 388 Array<Real,3>::RefArray<1> rQXY = rQX[m]; 389 #else 390 RefArray<Real,1> rGDXY = rGDX[j0]; 391 RefArray<Real,1> rVXY = rVX[m]; 392 RefArray<Real,1> rQXY = rQX[m]; 393 #endif 394 for (; n < nz; n++){ 395 int k0 = k-n; 396 if(pbcZ){ 397 k0 = (k0+nz)%nz; 398 k0 = std::min(k0,nz-k0); 399 } 400 else { 401 k0 = k0 < 0 ? -k0 : k0; 402 } 403 Real g = rGDXY[k0]; 404 v += rQXY[n]*g; 405 rVXY[n] += q*g; 406 #ifdef DEBUG_MULTIGRID_TIMING 407 myCounterDirect++; 408 #endif 409 } 410 } 411 } 412 rV1XY[k] += q*g0+v; 413 count++; 414 if(count >= size) 415 return; 416 } 417 } 418 } 419 } 420 421 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> coarseToFine(int level,unsigned int block,unsigned int nBlocks)422 inline void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::coarseToFine(int level, unsigned int block, unsigned int nBlocks){ 423 int nx = myDim[level-1].x; 424 int ny = myDim[level-1].y; 425 int nz = myDim[level-1].z; 426 int nx2 = myDim[level].x; 427 int ny2 = myDim[level].y; 428 int nz2 = myDim[level].z; 429 Array<Real,3>& rV0 = myV[level-1]; 430 Array<Real,3>& rV1 = myV[level]; 431 432 int nyz = ny*nz; 433 int n = nx*nyz; 434 int sn = (n*block)/nBlocks; 435 int en = (n*(block+1))/nBlocks - 1; 436 int count = 0; 437 int size = en-sn+1; 438 if(size == 0) 439 return; 440 441 int k = sn % nz; 442 int j = (sn / nz) % ny; 443 int i = (sn / nyz); 444 int ez = (en % nz)+1; 445 int ey = ((en / nz) % ny)+1; 446 int ex = (en / nyz)+1; 447 if(j < ey-1) 448 ez = nz; 449 if(i < ex-1){ 450 ey = ny; 451 ez = nz; 452 } 453 454 for (; i < ex; i++,j=0){ 455 int i1 = i/myRatio; 456 Real* thetaX = myGridInterpolation[i%myRatio].theta; 457 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 458 Array<Real,3>::RefArray<2> rV0X = rV0[i]; 459 #else 460 RefArray<Real,2> rV0X = rV0[i]; 461 #endif 462 for (; j < ey; j++,k=0){ 463 int j1 = j/myRatio; 464 Real* thetaY = myGridInterpolation[j%myRatio].theta; 465 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 466 Array<Real,3>::RefArray<1> rV0XY = rV0X[j]; 467 #else 468 RefArray<Real,1> rV0XY = rV0X[j]; 469 #endif 470 for (; k < ez; k++){ 471 int k1 = k/myRatio; 472 Real* thetaZ = myGridInterpolation[k%myRatio].theta; 473 Real v = 0.0; 474 for(int i0=0;i0<myInterOrder;i0++){ 475 int i2 = i1+i0; 476 if(pbcX) i2 = (myNXOffset+i2) % nx2; 477 Real x = thetaX[i0]; 478 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 479 Array<Real,3>::RefArray<2> rV1X = rV1[i2]; 480 #else 481 RefArray<Real,2> rV1X = rV1[i2]; 482 #endif 483 for(int j0=0;j0<myInterOrder;j0++){ 484 int j2 = j1+j0; 485 if(pbcY) j2 = (myNYOffset+j2) % ny2; 486 Real xy = x*thetaY[j0]; 487 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 488 Array<Real,3>::RefArray<1> rV1XY = rV1X[j2]; 489 #else 490 RefArray<Real,1> rV1XY = rV1X[j2]; 491 #endif 492 for(int k0=0;k0<myInterOrder;k0++){ 493 int k2 = k1+k0; 494 if(pbcZ) k2 = (myNZOffset+k2) % nz2; 495 v += rV1XY[k2]*xy*thetaZ[k0]; 496 } 497 } 498 } 499 rV0XY[k] += v; 500 count++; 501 if(count >= size) 502 return; 503 } 504 } 505 } 506 } 507 508 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> correction(int level,unsigned int block,unsigned int nBlocks)509 inline void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::correction(int level, unsigned int block, unsigned int nBlocks){ 510 Real scale = 1.0/myScale[level]; 511 int nx = myDim[level].x; 512 int ny = myDim[level].y; 513 int nz = myDim[level].z; 514 Array<Real,3>& rV = myV[level]; 515 Array<Real,3>& rQ = myQ[level]; 516 517 int nyz = ny*nz; 518 int n = nx*nyz; 519 int sn = 0; 520 int en = n - 1; 521 if(nBlocks > 1){ 522 if(myBlocksCorrection.empty() || 523 myBlocksCorrection[level].size() != nBlocks+1 || 524 myBlocksCorrection[level][myBlocksCorrection[level].size()-1] != en +1) 525 blocksCorrection(nBlocks); 526 sn=myBlocksCorrection[level][block]; 527 en=myBlocksCorrection[level][block+1]-1; 528 } 529 int count = 0; 530 int size = en-sn+1; 531 if(size == 0) 532 return; 533 534 int k = sn % nz; 535 int j = (sn / nz) % ny; 536 int i = (sn / nyz); 537 int ez = (en % nz)+1; 538 int ey = ((en / nz) % ny)+1; 539 int ex = (en / nyz)+1; 540 if(j < ey-1) 541 ez = nz; 542 if(i < ex-1){ 543 ey = ny; 544 ez = nz; 545 } 546 547 for (; i < ex; i++,j=0){ 548 int hi_l = i+myGCorrDimX; 549 int so_l = -myGCorrDimX+1; 550 int lo_l = i+so_l; 551 if(!pbcX){ 552 hi_l = std::min(hi_l,nx); 553 so_l+=-std::min(lo_l,0); 554 lo_l = std::max(lo_l,0); 555 } 556 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 557 Array<Real,3>::RefArray<2> rVX = rV[i]; 558 #else 559 RefArray<Real,2> rVX = rV[i]; 560 #endif 561 for (; j < ey; j++,k=0){ 562 int hi_m = j+myGCorrDimY; 563 int so_m = -myGCorrDimY+1; 564 int lo_m = j+so_m; 565 if(!pbcY){ 566 hi_m = std::min(hi_m,ny); 567 so_m+=-std::min(lo_m,0); 568 lo_m = std::max(lo_m,0); 569 } 570 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 571 Array<Real,3>::RefArray<1> rVXY = rVX[j]; 572 #else 573 RefArray<Real,1> rVXY = rVX[j]; 574 #endif 575 for (; k < ez; k++){ 576 int hi_n = k+myGCorrDimZ; 577 int so_n = -myGCorrDimZ+1; 578 int lo_n = k+so_n; 579 if(!pbcZ){ 580 hi_n = std::min(hi_n,nz); 581 so_n+=-std::min(lo_n,0); 582 lo_n = std::max(lo_n,0); 583 } 584 Real v = 0.0; 585 for (int l = lo_l, l2 = so_l; l < hi_l; l++,l2++){ 586 int l0 = l; 587 if(pbcX) l0 = (l0+nx)%nx; 588 int id = std::max(-l2,l2); 589 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 590 Array<Real,3>::RefArray<2> rQX = rQ[l0]; 591 Array<Real,3>::RefArray<2> rGCX = myGCorrection[id]; 592 #else 593 RefArray<Real,2> rQX = rQ[l0]; 594 RefArray<Real,2> rGCX = myGCorrection[id]; 595 #endif 596 for (int m = lo_m, m2 = so_m; m < hi_m; m++,m2++){ 597 int m0 = m; 598 if(pbcY) m0 = (m0+ny)%ny; 599 int jd = std::max(-m2,m2); 600 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 601 Array<Real,3>::RefArray<1> rQXY = rQX[m0]; 602 Array<Real,3>::RefArray<1> rGCXY = rGCX[jd]; 603 #else 604 RefArray<Real,1> rQXY = rQX[m0]; 605 RefArray<Real,1> rGCXY = rGCX[jd]; 606 #endif 607 for (int n = lo_n, n2 = so_n; n < hi_n; n++,n2++){ 608 int n0 = n; 609 if(pbcZ) n0 = (n0+nz)%nz; 610 int kd = std::max(-n2,n2); 611 v += rQXY[n0]*rGCXY[kd]; 612 #ifdef DEBUG_MULTIGRID_TIMING 613 myCounterCorrection++; 614 #endif 615 } 616 } 617 } 618 rVXY[k] += v*scale; 619 count++; 620 if(count >= size) 621 return; 622 } 623 } 624 } 625 } 626 627 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> energy(int level,unsigned int block,unsigned int nBlocks)628 inline Real MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::energy(int level, unsigned int block, unsigned int nBlocks){ 629 Real e = 0.0; 630 int nx = myDim[level].x; 631 int ny = myDim[level].y; 632 int nz = myDim[level].z; 633 634 int nyz = ny*nz; 635 int n = nx*nyz; 636 int sn = (n*block)/nBlocks; 637 int en = (n*(block+1))/nBlocks - 1; 638 639 int count = 0; 640 int size = en-sn+1; 641 if(size == 0) 642 return 0.0; 643 644 int k = sn % nz; 645 int j = (sn / nz) % ny; 646 int i = (sn / nyz); 647 int ez = (en % nz)+1; 648 int ey = ((en / nz) % ny)+1; 649 int ex = (en / nyz)+1; 650 if(j < ey-1) 651 ez = nz; 652 if(i < ex-1){ 653 ey = ny; 654 ez = nz; 655 } 656 657 Array<Real,3>& rV = myV[level]; 658 Array<Real,3>& rQ = myQ[level]; 659 for (; i < ex; i++,j=0){ 660 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 661 Array<Real,3>::RefArray<2> rVX = rV[i]; 662 Array<Real,3>::RefArray<2> rQX = rQ[i]; 663 #else 664 RefArray<Real,2> rVX = rV[i]; 665 RefArray<Real,2> rQX = rQ[i]; 666 #endif 667 for (; j < ey; j++,k=0){ 668 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 669 Array<Real,3>::RefArray<1> rVXY = rVX[j]; 670 Array<Real,3>::RefArray<1> rQXY = rQX[j]; 671 #else 672 RefArray<Real,1> rVXY = rVX[j]; 673 RefArray<Real,1> rQXY = rQX[j]; 674 #endif 675 for (; k < ez; k++){ 676 e += rVXY[k] * rQXY[k]; 677 count++; 678 if(count >= size) 679 return 0.5*e; 680 } 681 } 682 } 683 // Real *q = myQ[level].begin(); 684 // Real *v = myV[level].begin(); 685 // for(int i=sn;i<en+1;i++) 686 // e += v[i]*q[i]; 687 return 0.5*e; 688 } 689 690 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> sumQ(int level)691 inline Real MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::sumQ(int level){ 692 Real e = 0.0; 693 Array<Real,3>& rQ = myQ[level]; 694 for (int i = 0; i < myDim[level].x; i++){ 695 for (int j = 0 ; j < myDim[level].y; j++){ 696 for (int k = 0; k < myDim[level].z; k++){ 697 e += rQ[i][j][k]; 698 } 699 } 700 } 701 return e; 702 } 703 704 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> sumV(int level)705 inline Real MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::sumV(int level){ 706 Real e = 0.0; 707 Array<Real,3>& rV = myV[level]; 708 for (int i = 0; i < myDim[level].x; i++){ 709 for (int j = 0 ; j < myDim[level].y; j++){ 710 for (int k = 0; k < myDim[level].z; k++){ 711 e += rV[i][j][k]; 712 } 713 } 714 } 715 return e; 716 } 717 718 719 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> interpolateForce(Real q,unsigned int index,Vector3D & force)720 inline void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::interpolateForce(Real q, unsigned int index, Vector3D& force){ 721 int i0 = myScaledParticleIntPositions[index].x; 722 int j0 = myScaledParticleIntPositions[index].y; 723 int k0 = myScaledParticleIntPositions[index].z; 724 Real fx = 0.0; 725 Real fy = 0.0; 726 Real fz = 0.0; 727 Real* thetaX = myInterpolations[index].x.theta; 728 Real* thetaY = myInterpolations[index].y.theta; 729 Real* thetaZ = myInterpolations[index].z.theta; 730 Real* dThetaX = myInterpolations[index].x.dTheta; 731 Real* dThetaY = myInterpolations[index].y.dTheta; 732 Real* dThetaZ = myInterpolations[index].z.dTheta; 733 Array<Real,3>& rV = myV[0]; 734 for(int i=0;i<myInterOrder;i++){ 735 int i1 = i+i0; 736 if(pbcX) i1 = i1 % myNX; 737 Real ax = dThetaX[i]; 738 Real ay = thetaX[i]; 739 Real az = thetaX[i]; 740 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 741 Array<Real,3>::RefArray<2> rVX = rV[i1]; 742 #else 743 RefArray<Real,2> rVX = rV[i1]; 744 #endif 745 for(int j=0;j<myInterOrder;j++){ 746 int j1 = j+j0; 747 if(pbcY) j1 = j1 % myNY; 748 Real abx = ax * thetaY[j]; 749 Real aby = ay * dThetaY[j]; 750 Real abz = az * thetaY[j]; 751 #ifdef NO_PARTIAL_TEMPLATE_SPECIALIZATION 752 Array<Real,3>::RefArray<1> rVXY = rVX[j1]; 753 #else 754 RefArray<Real,1> rVXY = rVX[j1]; 755 #endif 756 for(int k=0;k<myInterOrder;k++){ 757 int k1 = k+k0; 758 if(pbcZ) k1 = k1 % myNZ; 759 Real term = rVXY[k1]; 760 fx -= term*abx*thetaZ[k]; 761 fy -= term*aby*thetaZ[k]; 762 fz -= term*abz*dThetaZ[k]; 763 } 764 } 765 } 766 force += Vector3D(fx*q*myHXr,fy*q*myHYr,fz*q*myHZr); 767 } 768 769 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> precomputeG()770 inline void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::precomputeG(){ 771 if(!myGDirect.resize(ArraySizes(myDim[myLevels-1].x)(myDim[myLevels-1].y)(myDim[myLevels-1].z))) 772 Report::report << Report::error 773 <<"[MultiGrid<>::precomputeG] Could not allocate memory for GDirect[" 774 775 <<myDim[myLevels-1].x<<"]["<<myDim[myLevels-1].y<<"]["<<myDim[myLevels-1].z<<"]." 776 <<Report::endr; 777 Real scale = 1.0/myScale[myLevels-1]; 778 for (int i = 0; i < myDim[myLevels-1].x; i++){ 779 for (int j = 0 ; j < myDim[myLevels-1].y; j++){ 780 for (int k = 0; k < myDim[myLevels-1].z; k++){ 781 Real r = sqrt(power<2>(i*myHX)+power<2>(j*myHY)+power<2>(k*myHZ)); 782 myGDirect[i][j][k] = TKernel::smoothKernel(r,myS,myRS)*scale; 783 } 784 } 785 } 786 787 // No correction if just one level specified. 788 if(myLevels < 2 ) 789 return; 790 myGCorrDimX = std::min((int)ceil(myRatio*myS/myHX),myNX); 791 myGCorrDimY = std::min((int)ceil(myRatio*myS/myHY),myNY); 792 myGCorrDimZ = std::min((int)ceil(myRatio*myS/myHZ),myNZ); 793 if(!myGCorrection.resize(ArraySizes(myGCorrDimX)(myGCorrDimY)(myGCorrDimZ))) 794 Report::report << Report::error <<"[MultiGrid<>::precomputeG] Could not allocate memory for GCorrection[" 795 <<myGCorrDimX<<"]["<<myGCorrDimY<<"]["<<myGCorrDimZ<<"]."<<Report::endr; 796 for (int i = 0; i < myGCorrDimX; i++){ 797 for (int j = 0 ; j < myGCorrDimY; j++){ 798 for (int k = 0; k < myGCorrDimZ; k++){ 799 Real r = sqrt(power<2>(i*myHX)+power<2>(j*myHY)+power<2>(k*myHZ)); 800 myGCorrection[i][j][k] = TKernel::smoothKernel(r,myS,myRS)-TKernel::smoothKernel(r,myS*myRatio,myRS/myRatio); 801 } 802 } 803 } 804 } 805 806 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> blocksCorrection(unsigned int nBlocks)807 inline void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::blocksCorrection(unsigned int nBlocks){ 808 809 if(myLevels < 2) 810 return; 811 myBlocksCorrection.resize(myLevels-1); 812 813 long nxyz = (myGCorrDimX*2-1)*(myGCorrDimY*2-1)*(myGCorrDimZ*2-1); 814 815 for(int level=0;level<myLevels-1;level++){ 816 int nx = myDim[level].x; 817 int ny = myDim[level].y; 818 int nz = myDim[level].z; 819 long size = 0; 820 for (int i = 0; i < nx; i++){ 821 int hi_l = i+myGCorrDimX; 822 int so_l = -myGCorrDimX+1; 823 int lo_l = i+so_l; 824 if(!pbcX){ 825 hi_l = std::min(hi_l,nx); 826 lo_l = std::max(lo_l,0); 827 } 828 for (int j = 0 ; j < ny; j++){ 829 int hi_m = j+myGCorrDimY; 830 int so_m = -myGCorrDimY+1; 831 int lo_m = j+so_m; 832 if(!pbcY){ 833 hi_m = std::min(hi_m,ny); 834 lo_m = std::max(lo_m,0); 835 } 836 for (int k = 0; k < nz; k++){ 837 int hi_n = k+myGCorrDimZ; 838 int so_n = -myGCorrDimZ+1; 839 int lo_n = k+so_n; 840 if(!pbcZ){ 841 hi_n = std::min(hi_n,nz); 842 lo_n = std::max(lo_n,0); 843 } 844 for (int l = lo_l; l < hi_l; l++){ 845 for (int m = lo_m; m < hi_m; m++){ 846 for (int n = lo_n; n < hi_n; n++){ 847 size++; 848 } 849 } 850 } 851 size -= nxyz; 852 } 853 } 854 } 855 //Report::report << allnodesserial<< plain << level<<":"<<(Real)nx*(Real)ny*(Real)nz*(Real)nxyz-(Real)size<<","<<size <<Report::endr; 856 857 std::vector<int> blocks; 858 blocks.push_back(0); 859 long count = 0; 860 long block = 0; 861 862 for (int i = 0; i < nx; i++){ 863 int hi_l = i+myGCorrDimX; 864 int so_l = -myGCorrDimX+1; 865 int lo_l = i+so_l; 866 if(!pbcX){ 867 hi_l = std::min(hi_l,nx); 868 lo_l = std::max(lo_l,0); 869 } 870 for (int j = 0 ; j < ny; j++){ 871 int hi_m = j+myGCorrDimY; 872 int so_m = -myGCorrDimY+1; 873 int lo_m = j+so_m; 874 if(!pbcY){ 875 hi_m = std::min(hi_m,ny); 876 lo_m = std::max(lo_m,0); 877 } 878 for (int k = 0; k < nz; k++){ 879 int hi_n = k+myGCorrDimZ; 880 int so_n = -myGCorrDimZ+1; 881 int lo_n = k+so_n; 882 if(!pbcZ){ 883 hi_n = std::min(hi_n,nz); 884 lo_n = std::max(lo_n,0); 885 } 886 for (int l = lo_l; l < hi_l; l++){ 887 for (int m = lo_m; m < hi_m; m++){ 888 for (int n = lo_n; n < hi_n; n++){ 889 count++; 890 } 891 } 892 } 893 count -= nxyz; 894 block++; 895 if((((Real)block*(Real)nxyz-(Real)count)*(Real)nBlocks)/((Real)nx*(Real)ny*(Real)nz*(Real)nxyz-(Real)size) >= (Real)blocks.size()) 896 blocks.push_back(block); 897 } 898 } 899 } 900 block = nx*ny*nz; 901 blocks[blocks.size()-1]=block; 902 903 //Report::report << plain << level<<": ("<<blocks.size() <<") "; 904 while(blocks.size()<= nBlocks) 905 blocks.push_back(block); 906 myBlocksCorrection[level] = blocks; 907 //Report::report <<myBlocksCorrection[level].size() <<" :"; 908 //for(int b=0;b<myBlocksCorrection[level].size();b++) 909 // Report::report << myBlocksCorrection[level][b]<<" "; 910 //Report::report <<Report::endr; 911 } 912 } 913 914 915 916 917 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> MultiGrid()918 MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::MultiGrid(): 919 myQ(0), 920 myV(0), 921 myGDirect(ArraySizes(0)(0)(0)), 922 myGCorrDimX(0), 923 myGCorrDimY(0), 924 myGCorrDimZ(0), 925 myMin(Vector3D(Constant::MAXREAL,Constant::MAXREAL,Constant::MAXREAL)), 926 myMax(Vector3D(-Constant::MAXREAL,-Constant::MAXREAL,-Constant::MAXREAL)), 927 myCMin(Vector3D(Constant::MAXREAL,Constant::MAXREAL,Constant::MAXREAL)), 928 myCMax(Vector3D(-Constant::MAXREAL,-Constant::MAXREAL,-Constant::MAXREAL)), 929 myH(Vector3D(0.0,0.0,0.0)), 930 myOrigin(Vector3D(0.0,0.0,0.0)), 931 myNX(0), 932 myNY(0), 933 myNZ(0), 934 myNXOffset(0), 935 myNYOffset(0), 936 myNZOffset(0), 937 myHX(0.0), 938 myHY(0.0), 939 myHZ(0.0), 940 myHXr(0.0), 941 myHYr(0.0), 942 myHZr(0.0), 943 myScaledParticleIntPositions(0), 944 myInterpolations(0), 945 myGridInterpolation(0), 946 myS(0.0), 947 myInterOrder(0), 948 myRatio(0), 949 myLevels(0), 950 mySigma(0) 951 #ifdef DEBUG_MULTIGRID_TIMING 952 , 953 myCounterDirect(0), 954 myCounterCorrection(0) 955 #endif 956 { 957 } 958 959 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> initialize(unsigned int n,Real s,int levels,int nx,int ny,int nz,int interOrder,int ratio,Vector3D min,Vector3D max,Vector3D h,Vector3D origin)960 void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::initialize(unsigned int n, Real s, int levels, 961 int nx, int ny, int nz, 962 int interOrder, int ratio, 963 Vector3D min, Vector3D max, 964 Vector3D h, Vector3D origin){ 965 966 //Report::report << Report::debug << s<<","<<levels <<","<< nx<<","<<ny <<","<<nz<<","<<interOrder<<","<<ratio<<","<<min<<","<<max<<","<<h<<","<<origin<<Report::endr; 967 968 //Report::report.reset(); 969 //Report::report.setf(std::ios::scientific); 970 //Report::report.precision(18); 971 // Frist do some fast checks on the input ... 972 if(s <= 0.0) 973 Report::report << Report::error << "[MultiGrid::initialize] s (="<<s<<") must be > 0."<<Report::endr; 974 975 if(levels < 1) 976 Report::report << Report::error << "[MultiGrid::initialize] levels (="<<levels<<") must be > 0."<<Report::endr; 977 978 if(ratio < 2) 979 Report::report << Report::error << "[MultiGrid::initialize] ratio (=" 980 <<ratio<<") must be > 1."<<Report::endr; 981 982 if(interOrder < 2 || interOrder % 2 != 0) 983 Report::report << Report::error << "[MultiGrid::initialize] order (=" 984 <<interOrder<<") must be > 1 and even."<<Report::endr; 985 986 if((pbcX && nx < interOrder) || (pbcY && ny < interOrder) || (pbcZ && nz < interOrder)){ 987 Report::report << Report::error << "[MultiGrid::initialize]"; 988 if(pbcX) Report::report << " nx (=" <<nx<<") must be > order."; 989 if(pbcY) Report::report << " ny (=" <<ny<<") must be > order."; 990 if(pbcZ) Report::report << " nz (=" <<nz<<") must be > order."; 991 Report::report<<Report::endr; 992 } 993 if((!pbcX && h.x <= 0.0) || (!pbcX && h.x <= 0.0) || (!pbcX && h.x <= 0.0)){ 994 Report::report << Report::error << "[MultiGrid::initialize]"; 995 if(pbcX) Report::report << " h.x (=" <<h.x<<") must be > 0."; 996 if(pbcY) Report::report << " h.y (=" <<h.y<<") must be > 0."; 997 if(pbcZ) Report::report << " h.h (=" <<h.z<<") must be > 0."; 998 Report::report<<Report::endr; 999 } 1000 1001 myScaledParticleIntPositions.resize(n); 1002 myInterpolations.resize(n,Interpolation3D(interOrder)); 1003 1004 myDim.resize(levels); 1005 myScale.resize(levels); 1006 myQ.resize(levels); 1007 myV.resize(levels); 1008 myBlocksCorrection.resize(0); 1009 1010 // Tells us if the interpolation is 1 for one theta and 0 for the others (when w=0), 1011 // such that we can save one interpolation point for the next grid. 1012 // We have to keep track of that for the allocation of arrays since we 1013 // will need an extra point such that we can keep the interpolation loops 1014 // easy. 1015 // Hermitian interpolation has this nice property. 1016 mySigma = (TInterpolation::isSigma(interOrder)?1:0); 1017 1018 myRatio = ratio; 1019 1020 // The interpolation between the grids. 1021 myGridInterpolation.resize(myRatio,TInterpolation(interOrder)); 1022 for(int i=0;i<myRatio;i++) 1023 myGridInterpolation[i].set((Real)i/(Real)myRatio); 1024 1025 myS = s; 1026 myRS = 1.0/s; 1027 myLevels = levels; 1028 myInterOrder = interOrder; 1029 myH = h; 1030 myOrigin = origin; 1031 1032 // Toplevel 1033 myDim[levels-1].x = nx; 1034 myDim[levels-1].y = ny; 1035 myDim[levels-1].z = nz; 1036 1037 1038 for(int i=0;i<levels;i++){ 1039 myScale[i] = power(ratio,i); 1040 } 1041 1042 1043 updateSize(min,max); 1044 1045 1046 Report::report << hint <<"MultiGrid: s="<<myS<<", ratio="<<myRatio<<", h="<<myHX 1047 <<","<<myHY<<","<<myHZ<<", "<<TInterpolation::keyword 1048 <<" "<<myInterOrder<<"th order "<<TKernel::keyword<<", levels="; 1049 for(int i=0;i<levels;i++) 1050 Report::report << "("<< myDim[i].x<<","<<myDim[i].y<<","<<myDim[i].z<<"),"; 1051 Report::report<<" finest grid="<<myMin<<"-"<<myMax<<", Gc("<<myGCorrDimX<<","<<myGCorrDimY<<","<<myGCorrDimZ<<")."<<Report::endr; 1052 1053 } 1054 1055 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> updateSize(Vector3D min,Vector3D max)1056 void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::updateSize(Vector3D min, Vector3D max){ 1057 1058 //Report::report.reset(); 1059 //Report::report.setf(std::ios::scientific); 1060 //Report::report.precision(4); 1061 Vector3D inMin(min),inMax(max); 1062 int nx=0; 1063 int ny=0; 1064 int nz=0; 1065 int nx1=0; 1066 int ny1=0; 1067 int nz1=0; 1068 int b = myInterOrder-mySigma-1; 1069 1070 // Vacuum: Minimal number of grid points of the finest grid 1071 if(!pbcX){ 1072 min.x -= myH.x*BORDER_TOLERANCE; 1073 max.x += myH.x*BORDER_TOLERANCE; 1074 nx1 = (int)(ceil((max.x-myOrigin.x)/myH.x)-floor((min.x-myOrigin.x)/myH.x)+ b+1); 1075 nx = nx1+BORDER; 1076 } 1077 if(!pbcY){ 1078 min.y -= myH.y*BORDER_TOLERANCE; 1079 max.y += myH.y*BORDER_TOLERANCE; 1080 ny1 = (int)(ceil((max.y-myOrigin.y)/myH.y)-floor((min.y-myOrigin.y)/myH.y) + b+1); 1081 ny = ny1+BORDER; 1082 } 1083 if(!pbcZ){ 1084 min.z -= myH.z*BORDER_TOLERANCE; 1085 max.z += myH.z*BORDER_TOLERANCE; 1086 nz1 = (int)(ceil((max.z-myOrigin.z)/myH.z)-floor((min.z-myOrigin.z)/myH.z)+ b+1); 1087 nz = nz1+BORDER; 1088 } 1089 // Vacuum: Compute number of grid points for the coarsest grid 1090 for(int i=1;i<myLevels;i++){ 1091 if(!pbcX) 1092 nx = (int)ceil((Real)(nx-1)/(Real)(myRatio))+ b + 1; 1093 if(!pbcY) 1094 ny = (int)ceil((Real)(ny-1)/(Real)(myRatio))+ b + 1; 1095 if(!pbcZ) 1096 nz = (int)ceil((Real)(nz-1)/(Real)(myRatio))+ b + 1; 1097 } 1098 1099 // Vacuum: No re-size if inside finest grid and 1100 // no major grid size changes 1101 if(!pbcX && !pbcY && !pbcZ){ 1102 if(min.x > myCMin.x && min.y > myCMin.y && min.z > myCMin.z && 1103 max.x < myCMax.x && max.y < myCMax.y && max.z < myCMax.z && 1104 nx+BORDER>=myDim[myLevels-1].x && 1105 ny+BORDER>=myDim[myLevels-1].y && 1106 nz+BORDER>=myDim[myLevels-1].z) 1107 return; 1108 } 1109 1110 if(pbcX) 1111 nx=myDim[myLevels-1].x; 1112 if(pbcY) 1113 ny=myDim[myLevels-1].y; 1114 if(pbcZ) 1115 nz=myDim[myLevels-1].z; 1116 1117 bool update = false; 1118 1119 if(pbcX || pbcY || pbcZ || 1120 nx1 > myNX || ny1 > myNY || nz1 > myNZ || 1121 fabs((float)(nx-myDim[myLevels-1].x))>BORDER || 1122 fabs((float)(ny-myDim[myLevels-1].y))>BORDER || 1123 fabs((float)(nz-myDim[myLevels-1].z))>BORDER){ 1124 1125 update = true; 1126 1127 // Toplevel, coarsest grid 1128 if(!myQ[myLevels-1].resize(ArraySizes(nx+(pbcX?0:mySigma))(ny+(pbcY?0:mySigma))(nz+(pbcZ?0:mySigma)))) 1129 Report::report << Report::error <<"[MultiGrid<>::initialize] Could not allocate memory for Q[" 1130 <<(nx+(pbcX?0:mySigma))<<"]["<<(ny+(pbcY?0:mySigma))<<"][" 1131 <<(nz+(pbcZ?0:mySigma))<<"]."<<Report::endr; 1132 if(!myV[myLevels-1].resize(ArraySizes(nx+(pbcX?0:mySigma))(ny+(pbcY?0:mySigma))(nz+(pbcZ?0:mySigma)))) 1133 Report::report << Report::error <<"[MultiGrid<>::initialize] Could not allocate memory for V[" 1134 <<(nx+(pbcX?0:mySigma))<<"]["<<(ny+(pbcY?0:mySigma))<<"][" 1135 <<(nz+(pbcZ?0:mySigma))<<"]."<<Report::endr; 1136 1137 myDim[myLevels-1].x = nx; 1138 myDim[myLevels-1].y = ny; 1139 myDim[myLevels-1].z = nz; 1140 1141 // The other levels 1142 for(int i=myLevels-2;i>=0;i--){ 1143 if(pbcX) 1144 nx = nx * myRatio; 1145 else 1146 nx = (nx-myInterOrder+mySigma) * myRatio + 1; 1147 1148 if(pbcY) 1149 ny = ny * myRatio; 1150 else 1151 ny = (ny-myInterOrder+mySigma) * myRatio + 1; 1152 1153 if(pbcZ) 1154 nz = nz * myRatio; 1155 else 1156 nz = (nz-myInterOrder+mySigma) * myRatio + 1; 1157 1158 if(!myQ[i].resize(ArraySizes(nx+(pbcX?0:mySigma))(ny+(pbcY?0:mySigma))(nz+(pbcZ?0:mySigma)))) 1159 Report::report << Report::error <<"[MultiGrid<>::initialize] Could not allocate memory for Q[" 1160 <<(nx+(pbcX?0:mySigma))<<"]["<<(ny+(pbcY?0:mySigma))<<"][" 1161 <<(nz+(pbcZ?0:mySigma))<<"]."<<Report::endr; 1162 if(!myV[i].resize(ArraySizes(nx+(pbcX?0:mySigma))(ny+(pbcY?0:mySigma))(nz+(pbcZ?0:mySigma)))) 1163 Report::report << Report::error <<"[MultiGrid<>::initialize] Could not allocate memory for V[" 1164 <<(nx+(pbcX?0:mySigma))<<"]["<<(ny+(pbcY?0:mySigma))<<"][" 1165 <<(nz+(pbcZ?0:mySigma))<<"]."<<Report::endr; 1166 myDim[i].x = nx; 1167 myDim[i].y = ny; 1168 myDim[i].z = nz; 1169 } 1170 1171 myNX = myDim[0].x; 1172 myNY = myDim[0].y; 1173 myNZ = myDim[0].z; 1174 1175 // Offset for mapping gridpoints between two levels 1176 myNXOffset = -(myInterOrder-1)/2; 1177 if(pbcX) myNXOffset += myNX; 1178 myNYOffset = -(myInterOrder-1)/2; 1179 if(pbcY) myNYOffset += myNY; 1180 myNZOffset = -(myInterOrder-1)/2; 1181 if(pbcZ) myNZOffset += myNZ; 1182 1183 // Mesh size of the finest grid 1184 if(pbcX) 1185 myHX = (max.x - min.x)/((Real)myNX); 1186 else 1187 myHX = myH.x; 1188 if(pbcY) 1189 myHY = (max.y - min.y)/((Real)myNY); 1190 else 1191 myHY = myH.y; 1192 if(pbcZ) 1193 myHZ = (max.z - min.z)/((Real)myNZ); 1194 else 1195 myHZ = myH.z; 1196 1197 myHXr = 1.0/myHX; 1198 myHYr = 1.0/myHY; 1199 myHZr = 1.0/myHZ; 1200 } 1201 1202 if(!pbcX){ 1203 min.x = floor((min.x-myOrigin.x)/myH.x-(myNX-nx1)/2)*myH.x+myOrigin.x; 1204 max.x = min.x+(myNX-1.0-b)*myH.x; 1205 myCMin.x = min.x; 1206 myCMax.x = max.x; 1207 min.x -= myH.x* (b/2); 1208 max.x += myH.x* (b-b/2); 1209 } 1210 if(!pbcY){ 1211 min.y = floor((min.y-myOrigin.y)/myH.y-(myNY-ny1)/2)*myH.y+myOrigin.y; 1212 max.y = min.y+(myNY-1.0-b)*myH.y; 1213 myCMin.y = min.y; 1214 myCMax.y = max.y; 1215 min.y -= myH.y* (b/2); 1216 max.y += myH.y* (b-b/2); 1217 } 1218 if(!pbcZ){ 1219 min.z = floor((min.z-myOrigin.z)/myH.z-(myNZ-nz1)/2)*myH.z+myOrigin.z; 1220 max.z = min.z+(myNZ-1.0-b)*myH.z; 1221 myCMin.z = min.z; 1222 myCMax.z = max.z; 1223 min.z -= myH.z* (b/2); 1224 max.z += myH.z* (b-b/2); 1225 } 1226 1227 if(!update && ((myMin-min).normSquared() > Constant::EPSILON || 1228 (myMax-max).normSquared() > Constant::EPSILON)){ 1229 Report::report.reset(); 1230 Report::report.precision(7); 1231 Report::report << Report::hint << "MultiGrid:"<<myMin<<"-"<<myMax<<" to "<<min<<"-"<<max<<Report::endr; 1232 } 1233 1234 myMin = min; 1235 myMax = max; 1236 1237 if(update){ 1238 myBlocksCorrection.resize(0); 1239 if(!myGDirect.empty()){ 1240 Report::report.reset(); 1241 Report::report.precision(7); 1242 Report::report << Report::hint << "MultiGrid: Finest grid("<<myNX<<","<<myNY<<","<<myNZ<<") was re-sized to "<<myMin<<"-"<<myMax<<" ("<<inMin<<"-"<<inMax<<")."<<Report::endr; 1243 } 1244 precomputeG(); 1245 } 1246 } 1247 1248 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> clear()1249 inline void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::clear(){ 1250 for(int level=0;level<myLevels;level++){ 1251 // Array<Real,3>& rV = myV[level]; 1252 // Array<Real,3>& rQ = myQ[level]; 1253 // for (int i = 0; i < myDim[level].x+(pbcX?0:mySigma); i++){ 1254 // for (int j = 0 ; j < myDim[level].y+(pbcY?0:mySigma); j++){ 1255 // for (int k = 0; k < myDim[level].z+(pbcZ?0:mySigma); k++){ 1256 // rV[i][j][k] = 0.0; 1257 // rQ[i][j][k] = 0.0; 1258 // } 1259 // } 1260 // } 1261 int n = myQ[level].size(); 1262 Real *q = myQ[level].begin(); 1263 Real *v = myV[level].begin(); 1264 for(int i=0;i<n;i++){ 1265 q[i] = 0.0; 1266 v[i] = 0.0; 1267 } 1268 } 1269 } 1270 1271 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> printQ(int level)1272 void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::printQ(int level){ 1273 Report::report << plain << "Level "<<level<<":"<<Report::endr; 1274 Real q = 0.0; 1275 1276 for (int i = 0; i < myDim[level].x; i++){ 1277 for (int j = 0 ; j < myDim[level].y; j++){ 1278 Report::report << plain << "Q["<<i<<"]["<<j<<"][0-"<<myDim[level].z-1<<"] : "; 1279 for (int k = 0; k < myDim[level].z; k++){ 1280 Report::report <<myQ[level][i][j][k]<<" "; 1281 q += myQ[level][i][j][k]; 1282 } 1283 Report::report <<Report::endr; 1284 } 1285 } 1286 Report::report << plain <<"Sum Q :"<<q<<Report::endr; 1287 } 1288 1289 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> printV(int level)1290 void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::printV(int level){ 1291 Report::report << plain << "Level "<<level<<":"<<Report::endr; 1292 Real v = 0.0; 1293 for (int i = 0; i < myDim[level].x; i++){ 1294 for (int j = 0 ; j < myDim[level].y; j++){ 1295 Report::report << plain << "V["<<i<<"]["<<j<<"][0-"<<myDim[level].z-1<<"] : "; 1296 for (int k = 0; k < myDim[level].z; k++){ 1297 Report::report <<myV[level][i][j][k]<<" "; 1298 v += myV[level][i][j][k]; 1299 } 1300 Report::report <<Report::endr; 1301 } 1302 } 1303 Report::report << plain <<"Sum V :"<<v<<Report::endr; 1304 } 1305 1306 template<class TInterpolation, class TKernel, bool pbcX, bool pbcY, bool pbcZ> printConst()1307 void MultiGrid<TInterpolation,TKernel,pbcX,pbcY,pbcZ>::printConst(){ 1308 for (int i = 0; i < myDim[myLevels-1].x; i++){ 1309 for (int j = 0 ; j < myDim[myLevels-1].y; j++){ 1310 Report::report << plain << "Gd["<<i<<"]["<<j<<"][0-"<<myDim[myLevels-1].z-1<<"] : "; 1311 for (int k = 0; k < myDim[myLevels-1].z; k++){ 1312 Report::report <<myGDirect[i][j][k]<<" "; 1313 } 1314 Report::report <<Report::endr; 1315 } 1316 } 1317 1318 // No correction if just one level specified. 1319 if(myLevels < 2 ) 1320 return; 1321 for (int i = 0; i < myGCorrDimX; i++){ 1322 for (int j = 0 ; j < myGCorrDimY; j++){ 1323 Report::report << plain << "Gc["<<i<<"]["<<j<<"][0-"<<myGCorrDimZ-1<<"] : "; 1324 for (int k = 0; k < myGCorrDimZ; k++){ 1325 Report::report <<myGCorrection[i][j][k]<<" "; 1326 } 1327 Report::report <<Report::endr; 1328 } 1329 } 1330 } 1331 } 1332 #endif 1333